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)
As per the philosophy of Monte Carlo and your attempt, I think you can try replicate + rowMeans (I think this combo applies to your update with fevd as well but I kept using optim here)
N <- 30
est <- replicate(N, {
x <- rgumbel(20, loc = 0, scale = 1)
m <- length(x)
fit <- optim(
par = c(1, 1),
fn = ll_fn, z = x, m = m,
hessian = TRUE, control = list(fnscale = -1)
)
fit$par
})
mse <- rowMeans((est - c(0, 1))^2)
and you will obtain
> est
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 0.03994695 0.0411953 0.04289633 0.4483689 -0.1757629 -0.002960991
[2,] 0.86737551 1.0367012 0.78748509 1.3409259 1.1527607 1.277664938
[,7] [,8] [,9] [,10] [,11] [,12]
[1,] -0.05394607 0.07246731 0.9500388 -0.2932722 -0.03107634 -0.01794733
[2,] 0.67182217 0.99549765 1.4256878 0.9107769 1.03009996 0.59352028
[,13] [,14] [,15] [,16] [,17] [,18] [,19]
[1,] -0.2268109 0.1411913 -0.1689498 -0.1452246 -0.3621658 0.02991549 0.1332464
[2,] 0.7631312 1.0881250 0.9639512 0.9912077 0.8454992 1.08480054 1.0938930
[,20] [,21] [,22] [,23] [,24] [,25] [,26]
[1,] 0.03264025 0.2279907 0.5108963 0.2505532 0.1572398 -0.4195196 -0.09234722
[2,] 0.89162837 1.0663945 1.0713617 0.8400802 0.6026962 0.8300022 0.61653388
[,27] [,28] [,29] [,30]
[1,] -0.1322638 0.0582058 -0.03012397 0.1494059
[2,] 1.0835487 0.8667572 1.05939105 0.8632393
> mse
[1] 0.07120414 0.04254979
In est, the first row is for mu, and the second row is for sigma, and mse consists of mu (1st val) and sigma (2nd val)