r代码运行始终报错Error in -t(X) %*% diag(pi_m) : 非整合参数，但是检查了行列数没有问题

Question

# 生成示例数据集
set.seed(0)
n <- 100  # 样本数量
p <- 3    # 自变量数量

X <- matrix(rnorm(n * p), ncol = p) #X为n*p

#intercept <- matrix(1, nrow = n, ncol = 1) #添加列为n*1
#X <- cbind(intercept, X)  # 添加截距列，X变为n*p+1


true_beta <- matrix(seq(1, p), ncol = 1)  # 生成从1到p的整数序列，并创建一个列向量p*1

pi <- 1 / (1 + exp(-X %*% true_beta))  # 计算概率，pi为n*1的矩阵
Y <- rbinom(n, 1, pi)  # 生成二进制结果，由于pi是n*1，所以Y为n*1的矩阵


# 牛顿-拉弗森方法的逻辑回归参数估计
newton_raphson <- function(X, Y, max_iter = 100, tol = 1e-6) {
  n <- nrow(X)  #X的行数
  p <- ncol(X)  #X的列数
  #beta_l <- rep(0, n)  # 初始化参数向量为0
  beta_l <- matrix(rep(0, p), ncol = 1)#p*1
  history <- list()  # 用于存储迭代历史记录
  
  

 
  for (i in 1:max_iter) {
    pi <- 1 / (1 + exp(-X %*% beta_l))  # 计算概率,pi为n*1
    pi_m <-  pi * (1 - pi)  
    gradient <- t(X) %*% (Y - pi)  # 计算一阶导数，t(X)为p*n，Y-pi为n*1
    hessian <- -t(X) %*% diag(pi_m) %*% X # 计算二阶导数，t(X)为p*n，diag()为n*n，X为n*p
  
    update <- solve(hessian) %*% gradient  # 计算参数更新量
    
    beta_new <- beta_l + update  # 更新参数向量
    history[[i]] <- beta_new  # 添加到历史记录
    
    if (sqrt(sum(update^2)) < tol) {  # 判断是否收敛
      break
    }
    
    beta_l <- beta_new
  }
  
  return(list(beta_l = beta_new, history = history))
}




# 调用函数进行参数估计
result <- newton_raphson(X, Y)
estimated_beta <- result$beta
history <- do.call(rbind, result$history)

# 绘制迭代历史记录
plot(history[, 1], type = "l", xlab = "Iterations", ylab = "Parameter Value",
     main = "Parameter Convergence", col = "red", ylim = range(history))
lines(history[, 2], col = "blue")
lines(history[, 3], col = "green")
lines(history[, 4], col = "purple")
legend("topright", legend = c("Intercept", "Beta 1", "Beta 2", "Beta 3"),
       col = c("red", "blue", "green", "purple"), lty = 1)

始终报错Error in -t(X) %*% diag(pi_m) : 非整合参数
但是反复检查了行列数，没有问题可以相乘呀
麻烦解释一下

Answer 1

【以下回答由 GPT 生成】

此错误通常是由于参数的维度不匹配引起的。在代码的hessian计算中，存在参数不匹配的问题。具体地说，在计算hessian时，它的维度应为pp，而不是像目前这样是np。

要解决这个问题，我们需要重新考虑计算hessian的方式。我们可以使用矩阵函数crossprod()来代替%*%运算符来计算t(X) %% diag(pi_m) %% X。这将确保我们得到正确维度的hessian矩阵。

下面是更新后的代码：

# 生成示例数据集
set.seed(0)
n <- 100  # 样本数量
p <- 3    # 自变量数量

X <- matrix(rnorm(n * p), ncol = p) #X为n*p

true_beta <- matrix(seq(1, p), ncol = 1)  # 生成从1到p的整数序列，并创建一个列向量p*1

pi <- 1 / (1 + exp(-X %*% true_beta))  # 计算概率，pi为n*1的矩阵
Y <- rbinom(n, 1, pi)  # 生成二进制结果，由于pi是n*1，所以Y为n*1的矩阵


# 牛顿-拉弗森方法的逻辑回归参数估计
newton_raphson <- function(X, Y, max_iter = 100, tol = 1e-6) {
  n <- nrow(X)  #X的行数
  p <- ncol(X)  #X的列数
  beta_l <- matrix(rep(0, p), ncol = 1)#p*1
  history <- list()  # 用于存储迭代历史记录

  for (i in 1:max_iter) {
    pi <- 1 / (1 + exp(-X %*% beta_l))  # 计算概率,pi为n*1
    pi_m <-  pi * (1 - pi)  
    gradient <- t(X) %*% (Y - pi)  # 计算一阶导数,t(X)为p*n，Y-pi为n*1
    hessian <- -crossprod(X, diag(pi_m), X) # 计算二阶导数,t(X)为p*n，diag()为n*n，X为n*p

    update <- solve(hessian) %*% gradient  # 计算参数更新量

    beta_new <- beta_l + update  # 更新参数向量
    history[[i]] <- beta_new  # 添加到历史记录

    if (sqrt(sum(update^2)) < tol) {  # 判断是否收敛
      break
    }

    beta_l <- beta_new
  }

  return(list(beta_l = beta_new, history = history))
}

# 调用函数进行参数估计
result <- newton_raphson(X, Y)
estimated_beta <- result$beta_l
history <- do.call(rbind, result$history)

# 绘制迭代历史记录
plot(history[, 1], type = 'l', xlab = 'Iterations', ylab = 'Parameter Value',
     main = 'Parameter Convergence', col = 'red', ylim = range(history))
lines(history[, 2], col = 'blue')
lines(history[, 3], col = 'green')
lines(history[, 4], col = 'purple')
legend('topright', legend = c('Intercept', 'Beta 1', 'Beta 2', 'Beta 3'),
       col = c('red', 'blue', 'green', 'purple'))

请尝试使用更新后的代码运行，应该不再出现参数维度不匹配的错误。

如果你已经解决了该问题, 非常希望你能够分享一下解决方案, 写成博客, 将相关链接放在评论区, 以帮助更多的人 ^-^