Search code examples
pythonnumpysimulationprobability

Why doesn't the following code give return a probability of approximately 0.75?


Tasked with simulating 10,000 throws of a pair of die and calculating how many of these throws produce an even number when we take the product of the results of the two dice. The idea is to show that this should be quite close to the theoretical probability of 0.75.

I wrote the following code but it gives me 8167 even-product throws when it should be something closer to 7500.

np.random.seed(193)
#np.random.randint(0,7) is a (random) die 

count=0
for i in range(10000):
    
if np.mod(np.random.randint(0,7)*np.random.randint(0,7), 2)==0: 
        count+=1

count

(I understand that there are many ways to do this and perhaps more elegant ones, just trying to see why this gave the result it did.)


Solution

  • As pointed out in the comments, you want np.random.randint(1, 7), as there's no 0 on dice:

    import numpy as np
    
    np.random.seed(193)
    
    count = 0
    for i in range(10000):
        if np.mod(np.random.randint(1, 7) * np.random.randint(1, 7), 2) == 0:
            count += 1
    
    print(count)
    

    Or just:

    import numpy as np
    
    np.random.seed(193)
    
    count = sum([1 - np.mod(np.random.randint(1, 7) * np.random.randint(1, 7), 2)
                 for _ in range(10000)])
    
    print(count)