STAN'da maksimum olabilirlik tahminleriyle ilgili standart hataları nasıl alabilirim?

Question

Ben Stan maksimum olabilirlik optimizasyonu kullanıyorum, fakat maalesef optimizing() fonksiyonu standart hataları bildirmez: Ben standart tahminleri hataları (veya güven aralığı - kantilleri) alırım nasıl

> 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 

ve muhtemelen p değerleri?

DÜZENLEME: Ben nazikçe @Ben Goodrich tarafından tavsiye etmedi:

> 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 

Fakat bu "kısıtsız" standart hatalar gerçek olanlardan çok farklı görünmektedir - burada Bayes uymasından çıkışı olarak stan() kullanarak: Sen çıkış listesinin bir parçası olarak Hessian dönecektir optimizing fonksiyonu için hessian = TRUE argüman belirtebilirsiniz

> 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 

Answer 1

. Böylece, tahmini standart hataları sqrt(diag(solve(-MLb4c$hessian))) aracılığıyla edinebilirsiniz; Ancak bu standart hatalar, sınırsız boşluğundaki tahminlere aittir. kısıtlı boşluğundaki parametreler için tahmini standart hataları elde etmek için, delta yöntemini kullanabilir veya ortalama vektörü MLb4c$par olan ve varyansı kovaryansı solve(-MLb4c$hessian) olan, çok fazla değişkenli normal dağılımdan çok kez çizebilir, constrain_pars işlevi ile boşluk ve her sütunun standart sapmasını tahmin edin. İşte

davanız

# 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)