Testing Monte Carlo Simulation
What is Monte Carlo Simulation?
Monte Carlo Simulation, also known as the Monte Carlo Method or a multiple probability simulation, is a mathematical technique, which is used to estimate the possible outcomes of an uncertain event. The Monte Carlo Method was invented by John von Neumann and Stanislaw Ulam during World War II to improve decision making under uncertain conditions.
Monte Carlo Simulations are utilized for long-term predictions due to their accuracy.
For eg:- Calculating probability of rolling 2 dice 10,000 times accurately.
The probability of an event \(A\) can be estimated as follows. We can simulate the experiment repeatedly and independently, say \(N\) times, and count the number of times the event occurred, say \(N_A\).
A good estimate of \(P(A)\) is the following:
$$P(A) \approx \frac{N_A}{N}$$
As \(N\) grows larger and larger, the estimate becomes better and better. This method is generally termed as Monte Carlo simulation.
Polya's urn scheme
Suppose an urn contains \(r\) red and \(b\) blue balls. The experiment proceeds in multiple steps, where Step \(i\) is as follows:
Step \(i\): Draw a ball at random, note down its colour and replace it in the urn. Add \(c\) more balls of the same colour to the urn.
Let \(R_i\) be the event that the \(i\)-th ball drawn is red. Let \(B_i\) be the event that the \(i\)-th ball drawn is black.
Clearly, \(P(R_1) = \frac{r}{r+b}\) and \(P(B_1)=\frac{b}{r+b}\). It is perhaps surprising that, irrespective of \(c\), we have, for all \(i\),
$$P(R_i) = \frac{r}{r+b}, P(B_i) = \frac{b}{r+b}.$$
To prove the above, you can use induction. Assume that the above is true for \(i\) and show it is true for \(i+1\). Starting with \(i=1\), by induction, the statement becomes true.
We will setup a Monte Carlo simulation for verifying \(P(R_i)\) above for a few steps.
What we have to prove?
We have to show original probability is equal to Monte Carlo Simulated probability.
Setup
* Import Library
import numpy as np
* Function for uniform outcome
\(n\) : number of outcomes in the sample space
Output : \(m\) outcomes selected uniformly at random from 1 to \(n\)
def uniform(n, m):
return np.random.randint(1, n+1, size = m)
* Simulation
\(N=\) Number of time experiment perform.
I perform this simulation for :
\(N=10,100,1000,100000,1000000\)
> N=10
# Repeating experiment 10 times i.e; N=10.
no = 0 #variable for storing number of event occurence
r = 10; b = 5 #assume 1 to r is red and r+1 to r+b is blue
print("Original Probability :",round(r/(r+b),3))
for i in range(10):
r = 10; b = 5
c = 3
for j in range(5): #do 5 steps
if uniform(r+b, 1) <= r:
r = r + c
else:
b = b + c
if uniform(r+b, 1) <= r: #in the 6th step, count if red ball drawn
no = no + 1
print("Simulated Probability:",no/10) #probability estimate by Monte Carlo
Output:-
Original Probability : 0.667
Monte Carlo Simulated Probability : 1.0
> N=100
# Repeating experiment 100 times i.e; N=100.
no = 0 #variable for storing number of event occurence
r = 10; b = 5 #assume 1 to r is red and r+1 to r+b is blue
print("Original Probability :",round(r/(r+b),3))
for i in range(100):
r = 10; b = 5
c = 3
for j in range(5): #do 5 steps
if uniform(r+b, 1) <= r:
r = r + c
else:
b = b + c
if uniform(r+b, 1) <= r: #in the 6th step, count if red ball drawn
no = no + 1
print("Monte Carlo Simulated Probability:",no/100) #probability estimate by Monte Carlo
Output:-
Original Probability : 0.667
Monte Carlo Simulated Probability : 0.73
> N=1000
# Repeating experiment 1000 times i.e; N=1000.
no = 0 #variable for storing number of event occurence
r = 10; b = 5 #assume 1 to r is red and r+1 to r+b is blue
print("Original Probability :",round(r/(r+b),3))
for i in range(1000):
r = 10; b = 5
c = 3
for j in range(5): #do 5 steps
if uniform(r+b, 1) <= r:
r = r + c
else:
b = b + c
if uniform(r+b, 1) <= r: #in the 6th step, count if red ball drawn
no = no + 1
print("Monte Carlo Simulated Probability:",no/1000) #probability estimate by Monte Carlo
Output:-
Original Probability : 0.667
Monte Carlo Simulated Probability : 0.664
> N=100000
# Repeating experiment 100000 times i.e; N=100000.
no = 0 #variable for storing number of event occurence
r = 10; b = 5 #assume 1 to r is red and r+1 to r+b is blue
print("Original Probability :",round(r/(r+b),3))
for i in range(100000):
r = 10; b = 5
c = 3
for j in range(5): #do 5 steps
if uniform(r+b, 1) <= r:
r = r + c
else:
b = b + c
if uniform(r+b, 1) <= r: #in the 6th step, count if red ball drawn
no = no + 1
print("Monte Carlo Simulated Probability:",no/100000) #probability estimate by Monte Carlo
Output:-
Original Probability : 0.667
Monte Carlo Simulated Probability : 0.66508
> N=1000000
# Repeating experiment 1000000 times i.e; N=1000000.
no = 0 #variable for storing number of event occurence
r = 10; b = 5 #assume 1 to r is red and r+1 to r+b is blue
print("Original Probability :",round(r/(r+b),3))
for i in range(1000000):
r = 10; b = 5
c = 3
for j in range(5): #do 5 steps
if uniform(r+b, 1) <= r:
r = r + c
else:
b = b + c
if uniform(r+b, 1) <= r: #in the 6th step, count if red ball drawn
no = no + 1
print("Monte Carlo Simulated Probability:",no/1000000) #probability estimate by Monte Carlo
Output:-
Original Probability : 0.667
Monte Carlo Simulated Probability : 0.66651
Conclusion
We see from above simulations that Original probability \(\approx\) Monte Carlo probability.
Hence, we can conclude that on performing an experiment like tossing coin, etc. for large number of times Monte Carlo estimate the possible outcomes of an uncertain event approximately correct.
Thank you for viewing my simulation.
\(Note:-\) I Use Python Programming Language for giving simulation.
