Optimizing function: incorrect estimation of the standard deviation of MLE

General
1

Hi,
I can’t get the correct standard deviation of MLEs with optimizing function in Stan. For example, We have the following formulas:

formula2

In following code, you can that optim function in R (stats package) works fine and optimizing function does not. Please guide me!

-Bistoon

library(rstan)
cat(
"
data{
int<lower=1> n;
vector[n] x;
}
parameters{
real<lower=0> sigma;
}
model{
x ~ normal(0, sigma);
}
"
, file = “model_Normal.stan”)

Stan_Normal = stan_model(“model_Normal.stan”)

generate data

set.seed(100)
n = 200
x = rnorm(n, 0, 2)

exact mles and sd of mles

exact_mle = sqrt(sum(x^2)/n)
exact_sd_mle = exact_mle * sqrt( 1 - (2/n)*(gamma((n+1)/2)/gamma(n/2))^2)

optim function in R (stats)

loglike = function(x, par) -sum(log(dnorm(x,0,par)))
lower = 0
upper = Inf
fit1 = optim(par= sd(x), loglike, x=x , hessian = TRUE, method = “L-BFGS-B”, lower=lower, upper=upper)
optim_mle= fit1$par;
optim_sd_mle = sqrt(diag(abs(solve(-(fit1$hessian)))))

optimizing function in rstan

fit2 = optimizing(Stan_Normal, hessian = TRUE, data=list(n=n, x=x), init = list(sigma = sd(x)))
optimizing_mle = fit2$par;
optimizing_sd_mle = sqrt(diag(abs(solve(-(fit2$hessian)))))

all the results

exact_mle
optim_mle
optimizing_mle

exact_sd_mle
optim_sd_mle
optimizing_sd_mle

2

Try optim with the log standard deviation as the parameter.

3

Thank you. I could get optimizing result with optim. How can I get correct results with optimizing function now?

4

The optimizing result does yield the correct result for the parameters in the unconstrained space, which is the only thing that matters to Stan. If you want the estimated variance-covariance matrix of the MLEs in the constrained space, then you either have to do the delta method calculations yourself or set the draws argument to some huge number and call cov on the theta_tilde element of the resulting list.

1 Like
5

I thought it would be easy! Can you give me an example of this?
Finally, we need the sd of estimators, like sampling function in Stan.

6

I did

cov(optimizing(..., draws = 10000)$theta_tilde)
2 Likes
7

Thank you very much.