# Optimizing function: incorrect estimation of the standard deviation of MLE

#1

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

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

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

#7

Thank you very much.