我如何在STAN中获得最大似然估计的标准误差？

Question

我如何在STAN中获得最大似然估计的标准误差？

10

我正在使用Stan中的最大似然优化，但不幸的是optimizing()函数没有报告标准误差：

> MLb4c <- optimizing(get_stanmodel(fitb4c), data = win.data, init = inits)  
STAN OPTIMIZATION COMMAND (LBFGS)
init = user
save_iterations = 1
init_alpha = 0.001
tol_obj = 1e-012
tol_grad = 1e-008
tol_param = 1e-008
tol_rel_obj = 10000
tol_rel_grad = 1e+007
history_size = 5
seed = 292156286
initial log joint probability = -4038.66
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes 
      13      -2772.49  9.21091e-005     0.0135987     0.07606      0.9845       15   
Optimization terminated normally: 
  Convergence detected: relative gradient magnitude is below tolerance
> t2 <- proc.time()
> print(t2 - t1)
   user  system elapsed 
   0.11    0.19    0.74 
> 
> MLb4c
$par
       psi      alpha       beta 
 0.9495000  0.4350983 -0.2016895 

$value
[1] -2772.489

> summary(MLb4c)
      Length Class  Mode   
par   3      -none- numeric
value 1      -none- numeric

我如何获取估计的标准误差（或置信区间-分位数），以及可能的p值？

编辑： 我按@Ben Goodrich的建议进行了操作：

> MLb4cH <- optimizing(get_stanmodel(fitb4c), data = win.data, init = inits, hessian = TRUE)

> sqrt(diag(solve(-MLb4cH$hessian)))
       psi      alpha       beta 
0.21138314 0.03251696 0.03270493

但是这些“无约束”的标准误估计看起来与真实的标准误估计非常不同 - 这里是使用stan()进行贝叶斯拟合的输出结果：

> print(outb4c, dig = 5)
Inference for Stan model: tmp_stan_model.
3 chains, each with iter=500; warmup=250; thin=1; 
post-warmup draws per chain=250, total post-warmup draws=750.

             mean se_mean      sd        2.5%         25%         50%         75%       97.5% n_eff    Rhat
alpha     0.43594 0.00127 0.03103     0.37426     0.41578     0.43592     0.45633     0.49915   594 1.00176
beta     -0.20262 0.00170 0.03167    -0.26640    -0.22290    -0.20242    -0.18290    -0.13501   345 1.00402
psi       0.94905 0.00047 0.01005     0.92821     0.94308     0.94991     0.95656     0.96632   448 1.00083
lp__  -2776.94451 0.06594 1.15674 -2780.07437 -2777.50643 -2776.67139 -2776.09064 -2775.61263   308 1.01220

- Tomas

频率学派的标准误差和贝叶斯标准差似乎很好地匹配了您的alpha和beta参数，只是不适用于psi。 - daknowles

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Ben Goodrich · Accepted Answer

您可以在optimizing函数中指定hessian = TRUE参数，这将返回Hessian作为输出列表的一部分。因此，您可以通过sqrt(diag(solve(-MLb4c$hessian)))获得估计的标准误差；然而，这些标准误差与非约束空间中的估计有关。要获取约束空间中参数的估计标准误差，您可以使用Delta方法或多次从多元正态分布中绘制，其平均向量为MLb4c$par，方差-协方差矩阵为solve(-MLb4c$hessian)，将这些绘制转换为受限空间使用constrain_pars函数，并估计每列的标准偏差。

以下是一些R代码，您可以根据自己的情况进行调整。

# 1: Compile and save a model (make sure to pass the data here)
model <- stan(file="model.stan", data=c("N","K","X","y"), chains = 0, iter = 0)

# 2: Fit that model
fit <- optimizing(object=get_stanmodel(model), as_vector = FALSE,
                   data=c("N","K","X","y"), hessian = TRUE)

# 3: Extract the vector theta_hat and the Hessian for the unconstrained parameters
theta_hat <- unlist(fit$par)
upars <- unconstrain_pars(linear, relist(theta_hat, fit$par))
Hessian <- fit$hessian

# 4: Extract the Cholesky decomposition of the (negative) Hessian and invert
R <- chol(-Hessian)
V <- chol2inv(R)
rownames(V) <- colnames(V) <- colnames(Hessian)

# 5: Produce a matrix with some specified number of simulation draws from a multinormal
SIMS <- 1000
len <- length(theta_hat)
unconstrained <- upars + t(chol(V)) %*% 
  matrix(rnorm(SIMS * len), nrow = len, ncol = SIMS)
theta_sims <- t(apply(unconstrained, 2, FUN = function(upars) {
  unlist(constrain_pars(linear, upars))
}))

# 6: Produce estimated standard errors for the constrained parameters
se <- apply(theta_sims, 2, sd)