Tuesday, January 31, 2012

Monte Carlo Simulation in Python

Problem:
Suppose you flip coins for a certain times, say n. You want to know your chance of getting certain number of heads in a row, say k heads in a row.

Let's start with the function that'll do a certain number of coin flips:
Let's do 10 flips right away:
do_n_flips(10)
'0100111001'

Let's suppose you get 1 euro if you get the given number of heads in a row while doing a certain number of flips. Here's the function that'll calculate your payoff:
Let's try our payoff function with doing 10 flips and getting payed for 3 heads in a row. sometimes it is 0, sometimes it is 1:
payoff(10, 3)
0.0
payoff(10, 3)
1.0

Now, what is your chance to get your 1 euro for 3 heads in a row from 10 flips? This brings us to the Monte Carlo simulation, which I'wont describe here at all, just give you a function which answers our question regarding the chance:
Running with a million iteration gives us a good approximation for the chance:
monte_carlo_solve(10, 3, 1000000)
0.507753

This post was inspired by a similar post by Remis.

3 comments:

Ian Bicking said...

Instead of "do_n_flips(n).find(''.rjust(k, '0')) > -1" you should do: "('0' * k) in do_n_flips(n)"

Sakti Dwi Cahyono said...

here my iterations result:

1 = 1.0
10 = 0.7
100 = 0.51
1000 = 0.495
10000 = 0.5097
100000 = 0.50985
1000000 = 0.508069
10000000 = 0.5073828

Jon said...

Nice. The analytical answer (from a Markov chain method) is 0.5078125 (=520/1024).