Monte Carlo Method

Updated at 2017-07-11 13:06

Monte Carlo method: Pouring out a box of coins on a table, and then computing the ratio of coins that land heads versus tails.

Monte Carlo method is the brute force of probability mathematics. Monte Carlo method utilizes randomness to solve probabilistic problems. It's most useful when it is difficult or impossible to create probability formula. But can be used for simplicity or if you just feel lazy.

var heads = 0;
var tails = 0;
// 3 million trials should be enough, right?
for (var i = 0; i < 3000000; i++) {
    if (Math.random() <= 0.5) {
        heads++
    } else {
        tails++
    }
}
console.log(heads, tails);

Monte Carlo relies on the law of large numbers. With a decent random generator and enough trials, you can get probability of anything you can simulate.

Monte Carlo methods have five steps:

Define the range of possible inputs.
Generate inputs randomly; they might have different probabilities.
Perform deterministic computation or "check".
Loop steps 2 and 3 for a large number of trials.
Aggregate the results to get probabilities.

Example; Monte Carlo method vs. actual math. We want to get probability of drawing one of two specific cards in a 30 card deck where starting hand is 3 cards.

from copy import copy
from random import shuffle

trials = 30000
success_count = 0
cards_in_deck = list(range(0, 30))

for i in range(0, trials):
    deck = copy(cards_in_deck)
    shuffle(deck)
    hand = [deck.pop(), deck.pop(), deck.pop()]
    if any(card == 0 for card in hand) or any(card == 1 for card in hand):
        success_count += 1

print('{}%'.format(success_count / trials * 100))
# => around 19%, which is close what hypergeometric distribution gives us

from operator import mul
from functools import reduce
from fractions import Fraction

def combinations(n, k):
    return int(reduce(mul, (Fraction(n - i, i + 1) for i in range(k)), 1))

def hypergeometric(population, wins, sample_size, required_successes):
    losses = population - wins
    sample_successes = combinations(wins, required_successes)
    sample_losses = combinations(losses, sample_size - required_successes)
    all_combinations = combinations(population, sample_size)
    return sample_successes * sample_losses / all_combinations

print('{}%'.format((hypergeometric(30, 2, 3, 1) * 100)))
# => 18.620689655172416%
# which you can double check http://stattrek.com/online-calculator/hypergeometric.aspx
# => 18.6206896551724%

Monte Carlo Method

Sources