How should I compute the Monte Carlo mean square error in R?
I try to get the Monte Carlo mean square error of the maximal likelihood estimators in R. I can write the calculation for the MLE that is repeated once, but I need to repeat the Monte Carlo calculation many times. How should I write this in R?
Firstly, we have sample size n=20 and the repetition N=30. The idea is that for the first repetition, we sample n=20 data from the Gumbel distribution with location 0 and scale 1 (true value mu=0 and sigma=1).
library(VGAM)
set.seed(101)
x <- rgumbel(20, loc = 0, scale = 1)
Then we compute the MLE of two parameters mu and sigma by the likelihood function
# Log-likelihood function of Gumbel distribution
ll_fn <- function(par, z, m){
mu <- par[1]
sigma <- par[2]
-m * log(sigma) - sum((z - mu)/sigma) - sum(exp(-((z - mu)/sigma)))
}
m <- length(x)
fit<- optim(par = c(1, 1),
fn = ll_fn, z = y, m = m,
hessian = TRUE, control = list(fnscale = -1))
fit3$par
#[1] 0.1718418 0.9813649
The Monte Carlo mean square error is that for the MLE of each i-th repetition

Comments
Post a Comment