如何从R的optim()函数提取Hessian矩阵计算二分类模型观测信息矩阵

你当前代码的核心问题是optim函数的参数位置错误：hessian=TRUE被错误放在了optim的调用括号外，属于语法错误，无法正常返回Hessian矩阵。

修正步骤 #

• 第一步：调整optim的参数位置，把hessian=TRUE移入optim的参数列表，开启Hessian矩阵计算
• 第二步：从optim返回结果中提取Hessian矩阵，该矩阵是你定义的负对数后验函数在最优参数处的二阶导数矩阵
• 第三步：根据你给出的规则观测信息矩阵 = - Hessian，直接计算得到观测信息矩阵即可

修正后的完整函数示例 #

posterior.mode <- function(y, x)
{
  # 样本量
  n <- length(y)
  # beta参数数量
  d <- ncol(x)
  # 第i个样本对应的以beta向量为输入的对数后验函数
  log.like_one <- function(i, beta)
  {
    # 给定X_i和beta时的p_i值
    Phi.xb <- pnorm( sum(x[i,] * beta) )
    return( log(1 - Phi.xb) + y[i]*log( Phi.xb/(1 - Phi.xb) ) )
  }
  # 全样本对数似然
  loglike <- function(beta)
  {
    L <- 0
    for(ii in 1:n) L <- L + log.like_one(ii,beta)
    return(L)
  }
  # 全样本负对数后验
  log.posterior <- function(beta) -loglike(beta) + ((beta/2)-log(1+exp(beta)))
  # 寻找对数后验最优的beta值，初始值默认全为0，开启Hessian计算
  optim_res <- optim(rep(0,d), log.posterior, hessian = TRUE)
  # 提取Hessian
  hessian_mat <- optim_res$hessian
  # 计算观测信息矩阵
  obs_info_mat <- - hessian_mat
  # 可按需返回多个结果，比如同时返回最优参数、Hessian、观测信息矩阵
  return(list(
    opt_beta = optim_res$par,
    hessian = hessian_mat,
    obs_info = obs_info_mat
  ))
}

调用示例 #

# 假设你已经准备好输入数据y和x
res <- posterior.mode(y, x)
# 直接提取观测信息矩阵
my_obs_info <- res$obs_info

注意事项 #

你当前定义的log.posterior是负对数后验，optim返回的Hessian是负对数后验的二阶导数，如果你后续计算参数标准误等指标，注意确认公式和矩阵的对应关系，避免符号错误。

内容的提问来源于stack exchange，提问作者M.Salim