CODEWITHSUNDEEP
pythonmathmontecarlo
Laura Marinescu
Laura Marinescu

Reputation: 57

Monte Carlo simulation dices in python

I have a question related to Monte Carlo simulation based on the probability of rolling 2 dices. How can you present when coding in python the fact that the sum is larger than n and smaller than m? As an example I made this in mathlab:

NT = 10^5; %number of throws
log = zeros(1,12);
for throw = 1:NT
    dices = ceil(6*rand(1,2));
    s = sum(dices);
    log(s) = log(s)+1;
end
p = 36*log(6:9)/NT;
s1 = sum(round(p))

In the above example I presumed that n is 5 and m is 10.

Thank you

Upvotes: 1

Views: 3014

Answers (2)

Campbell McDiarmid
Campbell McDiarmid

Reputation: 494

In each loop you want to simulate two separate random dice throws. My following snippet uses a list (you can use a dict if you like) to store the results of NT simulations:

import random

num_throws = 10**5  # NT
roll_log = [0] * 12  # Generate list for dice roll tallies

for i in range(num_throws):
    # Random integer between 1 and 6 inclusive for each dice
    dice_1 = random.randint(1, 6)
    dice_2 = random.randint(1, 6)

    # Sum the random dice and increment the tally for that particular roll total
    roll_sum = dice_1 + dice_2
    roll_log[roll_sum-1] += 1  # minus 1 because Python is 0-indexed

To process your resulting data you can access the tally for a particular dice roll in the results list by roll_log[roll-1] where 2 <= roll <= 12 (roll = 1 will have a probability of zero since it is impossible with 2+ dice). The following for-loop is just an example of how you can access the results of NT simulations if you are unfamiliar with enumeration in Python:

for i, tally in enumerate(roll_log):
    roll_prob = float(tally) / num_throws  # Experimental probability of roll
    roll = i + 1  # Since Python lists are 0-indexed
    print('{}: {}/{} = {}'.format(roll, tally, num_throws, roll_prob))

Output:

1: 0 / 100000 = 0
2: 2741 / 100000 = 0.02741
3: 5518 / 100000 = 0.05518
4: 8202 / 100000 = 0.08202
5: 11235 / 100000 = 0.11235
6: 14046 / 100000 = 0.14046
7: 16520 / 100000 = 0.1652
8: 13799 / 100000 = 0.13799
9: 11025 / 100000 = 0.11025
10: 8459 / 100000 = 0.08459
11: 5672 / 100000 = 0.05672
12: 2783 / 100000 = 0.02783

Specifically addressing the last part of your question, finding the probability of the dice roll being between n = 5 and m = 10 not inclusive, this can be done using a method called list slicing:

n = 5
m = 10
                                 #    6      7      8      9
rolls_between = roll_log[n:m-1]  # [14046, 16520, 13799, 11025]
sum_rolls_between = sum(rolls_between)  # 55390
prob_between = float(sum_rolls_between) / num_throws  # 0.5539

Note: The float conversion of either sum_rolls_between or num_throws on the last line is essential to obtain a decimal output, because division between two integers in Python always results in an integer output after applying mathematical floor() function. In other words, without changing one of those two values to a floating point value, the result will be 0.

Upvotes: 1

Aritesh
Aritesh

Reputation: 2103

See below- 

import numpy as np
NT = 10**5
n=5
m=10
x = np.random.randint(1, 12, NT)
s = sum((x>=n) & (x<=m))
p = s*1.0/NT
print(p)

Upvotes: 1

Related Questions