While the last post dealt with Bootstrap Sampling, no sampling discussion can be complete without discussion ‘Monte Carlo’ simulation. Readers please note, I will *not **discuss “MCMC (Markov Chain Monte Carlo)” *(perhaps in the future). MCMC primarily deals with distribution of equilibrium of the given Markov Chain.

Monte Carlo simulation is much used (perhaps overused) sampling for all physical simulation techniques. While the algorithm for generating Monte Carlo samples is very simple, its usage is widespread. Particularly in Computational Physics for modeling behavior of sub-atomic particles, in Robotics for generating possible states where the actual game state could be exponentially large and in bioinformatics for probabilistic graphical models and last but not least Monte Carlo Integration.

Readers not familiar with Monte Carlo may benefit from this tutorial by Jonathan Pengelly, as this blogpost restricts to presenting the simple code in R to do Monte Carlo simulation. Having a clean and grounds-up code is always beneficial as this helps tweak and reformulate the basics.

Following is the output obtained from the code given below:

Output

Code

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43

#Monte Carlo simulation 
#aaditya - 26/10/2012
 
######################################
StartEq = 100
SuccessPcnt = 61.54
SuccessAvg = 2444.78
FailureAvg = 1667.37
Iterations = 1000
NoOfEvents = 100
Rf = 1
#####################################
eq<-rep(NA, NoOfEvents)
 
doMCRun <- function() {
  eq[1] <- StartEq
  for(i in 2:NoOfEvents) {
    if(runif(1,,100)< SuccessPcnt)
      pl= SuccessAvg * Rf
    else
      pl= -1*FailureAvg * Rf
      eq[i] = eq[i-1]  + pl
  }
  return(eq)
}
 
values <- replicate(Iterations, doMCRun())
 
par(xaxs="i")
plot(1:NoOfEvents, rep(NA, NoOfEvents ), 
     xlab="Iterations", ylab="Growth",
     ylim=c(min(values),max(values)),
     xlim=c(1,NoOfEvents), main="Monte Carlo Simulation")
matlines(values,type="l",lty=1)
 
sd1 = sd(values[NoOfEvents,])
mean1 = mean(values[NoOfEvents,])
print(summary(values[NoOfEvents,]))
sdString = paste("SD : ", sd1)
write(sdString, file="")
outputString = paste("With 99.73% confidence the final growth will be between", 
                     mean1 - 3*sd1, " and ", mean1 + 3*sd1) 
write(outputString, file="")