Binomial Distribution problem using Montee Carlo

Kakashi Sensei 0
Viewed 42

Here, I tried to compute probability of getting one 4 times when I toss a die 6 times. Here is my code,

no=0

import random as rd
for i in range(1000):
  l=[rd.randint(1,6) for i in range(6)]
  a=l.count(1)
  if a==4:
    no+=1

print(no/1000)

I want to know, If this is actually a correct Montee Carlo approach of Binomial problem? And Is it correct?

probability-distribution probability-density montecarlo probability statistics
1 Answers

Your code is correct! Let me make two remarks to illustrate in more details how you can verify yourself that the code is correct.

  1. You already mention that the probability is related to the binomial distribution. That is, if getting a 1 is considered a success, then you need exactly 4 successes in n=6 trials when the success probability p=1/6. This probability can be computed exactly using for instance binom in scipy.stats.
  2. The outcome of a Monte Carlo experiment will always depend on the random variables you draw. Theoretically, we can only recover the true probability if we could have an infinite number of Monte Carlo replications. This is of course impossible in practice. However, it is usually a good idea to explicitly define the number of replications in the code (e.g. I took MonteCarloTrials = 1000000 in the code below). This allows you to increase the number of Monte Carlo experiments as you desire.
import random as rd
from scipy.stats import binom   # import binomial distribution

# Monte Carlo computation
no = 0
MonteCarloTrials = 1000000
for i in range(MonteCarloTrials):
  l = [rd.randint(1,6) for i in range(6)]
  a = l.count(1)
  if a==4:
    no+=1

# Exact binomial computation
n, p = 6, 1/6
x = binom.pmf(4, n, p) # in n=6 trials with should 4 successes with succes probability 1/6


print('Approximate probability computed with Monte Carlo:', no/MonteCarloTrials)
print('Exact probability:', x)
Hanno Reuvers0
Related