Using the Improved Simulation Method by asintunado.wordpress.com
set.seed(9999)
N <- 1377
M <- 6
SIM <- 1000000
se <- sqrt(0.5 * 0.5 /N)
duterte <- rnorm(SIM, 0.27, se)
poe <- rnorm(SIM, 0.23, se)
binay <- rnorm(SIM, 0.20, se)
roxas <- rnorm(SIM, 0.18, se)
santiago <- rnorm(SIM, 0.03, se)
santiago[santiago < 0 ] <- 0
ewan <- rnorm(SIM, 0.09, se)
ewan[ewan <0] <- 0
simulation <- matrix(nrow = M, ncol = SIM)
for(i in 1:SIM){
simulation[, i] <- rmultinom(1, N, c(duterte[i],
poe[i],
binay[i],
roxas[i],
santiago[i],
ewan[i]))
}
e <- table(apply(simulation, MARGIN=2, FUN = which.max))
print(e/sum(e))
##
## 1 2 3 4
## 0.930828 0.066467 0.002565 0.000140