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