In mcomas/normalmultinomial: Normal-Multinomial Distribution

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

• Mean:

knitr::kable(t(sapply(res, function(r) setNames(apply(r,1,mean)[2:3], c('M1','M2')))), digits = 4)

• Standard deviation:

knitr::kable(t(sapply(res, function(r) setNames(apply(r,1,sd)[2:3], c('M1','M2')))), digits = 4)

mcomas/normalmultinomial documentation built on May 22, 2019, 3:15 p.m.