It was my understanding that MLE estimates were asymptotically unbiased. The following simulation therefore confuses me. Can anyone help me with my understanding here?
I estimate the parameters of 50 random Weibulls that I've simulated with shape = scale = 1 (using the Wikipedia/R parameterisation)
After lots of simulations, we should* get a distribution of parameter estimates, the median of which equals the simulated values.
The result is that the median of the estimates are significantly different from the simulated parameters:
shape scale
est 1.018636 0.9954484
2.5% 1.015124 0.9919716
97.5% 1.021469 0.9984734
I expect I may be wrong at the "should*" above. If so, what is correct?
The R code is below:
I am simulating using rweibull()
and estimating parameters using mle()
in the core stats4 package and also fitdist()
in the fitdistrplus. I am using the correct estimates as the starting parameters, so there shouldn't be any optimisation issues here. (?)
# fitdistrplus::fitdist() fitting
res <- replicate(1e4,
fitdist(
rweibull(50,shape=1,scale=1),
distr="weibull",
start=list(shape=1,scale=1)
)$est
)
# stats4::mle() fitting
res <- replicate(1e4,{
y <- rweibull(50,shape=1,scale=1)
nll <- function(shape,scale) -sum(dweibull(y,shape,scale,log=TRUE))
res <- stats4::mle(
nll,
start=list(shape=1,scale=1)
)@coef
})
# calculate median of the simulated estimates and its bootstrapped confidence interval
apply(res,1,function(x){
c(est=median(x),
quantile(
replicate(1e3,
median(sample(x,replace=TRUE))
),
c(0.025,0.975)
)
)
})