Hi,

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

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