library(normalmultinomial) library(DT) set.seed(1)
Es seleccionen els paràmetres:
x = c(1,0) MU = 0 SIGMA = matrix(1)
Es calcula el valor esperat amb cada mètode utilitzant 1000 simulacions.
NSIM = 1000 fit0 <- expectedMonteCarlo(x, MU, SIGMA, nsim = NSIM, antithetic_variates = F, importance_sampling_mu = F) fit1 <- expectedMonteCarlo(x, MU, SIGMA, nsim = NSIM, antithetic_variates = T, importance_sampling_mu = F) fit2 <- expectedMonteCarlo(x, MU, SIGMA, nsim = NSIM, antithetic_variates = F, importance_sampling_mu = T) fit3 <- expectedMonteCarlo(x, MU, SIGMA, nsim = NSIM, antithetic_variates = T, importance_sampling_mu = T) fit4 <- expectedMetropolis(x, MU, SIGMA, nsim = NSIM) par(mfrow=c(1,2)) a0 = fit0[[2]] a1 = fit1[[2]] a2 = fit2[[2]] a3 = fit3[[2]] a4 = fit4[[2]] plot(cumsum(a0)/seq_along(a0), type='l', xlab='Iterations', ylab = 'Estimation', main='First moment') lines(cumsum(a1)/seq_along(a1), col='blue') lines(cumsum(a2)/seq_along(a2), col='green') lines(cumsum(a3)/seq_along(a3), col='orange') lines(cumsum(a4)/seq_along(a4), col='red') legend('bottomright', c('Montecarlo', 'M_AV', 'M_IS', 'M_AV_IS','Metropolis'), col=c('black', 'blue', 'green', 'orange','red'), bty = 'n', lty=1) a0 = c(fit0[[3]]) a1 = c(fit1[[3]]) a2 = c(fit2[[3]]) a3 = c(fit3[[3]]) a4 = c(fit4[[3]]) plot(cumsum(a0)/seq_along(a0), type='l', xlab='Iterations', ylab = 'Estimation', main='Second moment') lines(cumsum(a1)/seq_along(a1), col='blue') lines(cumsum(a2)/seq_along(a2), col='green') lines(cumsum(a3)/seq_along(a3), col='orange') lines(cumsum(a4)/seq_along(a4), col='red') legend('bottomright', c('Montecarlo', 'M_AV', 'M_IS', 'M_AV_IS','Metropolis'), col=c('black', 'blue', 'green', 'orange','red'), bty = 'n', lty=1)
Es repeteix el proces 500 vegades amb cada mètode per mesurar-ne la variabilitat.
res = list() res[['Montecarlo']] = replicate(500, sapply( expectedMonteCarlo(x, MU, SIGMA, nsim = NSIM, antithetic_variates = F, importance_sampling_mu = F), mean)) res[['M_AV']] = replicate(500, sapply( expectedMonteCarlo(x, MU, SIGMA, nsim = NSIM, antithetic_variates = T, importance_sampling_mu = F), mean)) res[['M_IS']] = replicate(500, sapply( expectedMonteCarlo(x, MU, SIGMA, nsim = NSIM, antithetic_variates = F, importance_sampling_mu = T), mean)) res[['M_AV_IS']] = replicate(500, sapply( expectedMonteCarlo(x, MU, SIGMA, nsim = NSIM, antithetic_variates = T, importance_sampling_mu = T), mean)) res[['Metropolis']] = replicate(500, sapply( expectedMetropolis(x, MU, SIGMA, nsim = NSIM), mean))
knitr::kable(t(sapply(res, function(r) setNames(apply(r,1,mean)[2:3], c('M1','M2')))), digits = 4)
knitr::kable(t(sapply(res, function(r) setNames(apply(r,1,sd)[2:3], c('M1','M2')))), digits = 4)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.