CODEWITHSUNDEEP
pythonnumpymatplotlibsimulationscientific-computing
user4933
user4933

Reputation: 1585

Fair coin in Python

I am trying to simulate a fair coin flip as independent bernoulli trials where the variable either goes to 1 or 0 then take the average at each trial for separate sequences.

My problem is I am getting one sequence which looks wrong and diverges from the others (Seq_1). Can you see what I am doing wrong.

# Parameters

nSeq = 5
nTrials = 10**3
p=0.5 # Fair coin

# Containers

x = np.zeros(nTrials+1, float)
nrange = range(nTrials+1)

# Simulation

for j in range(nSeq): 
    Mean_list = [0]  # re-initialize mean list for each sequence
    for i in range(nTrials):
        x[i+1] =  np.random.binomial(1, p) 
        xbar = np.mean(x)
        Mean_list.append(xbar) 
        
    plt.plot(nrange,Mean_list,label='Seq_'+str(j+1))
plt.ylabel('Estimate of p')
plt.xlabel('Trials')
plt.legend(loc=0,ncol=3,fontsize='small')
plt.show()

Which then gave me this code which does not look correct. It is correct to converge to 0.5 as the coin is fair but why does seq_1 always come out different to the others?

enter image description here

**This is trying to reproduce a graph from a book below

enter image description here

Summary

So to summarize. How does one generate a fair coin flip in python (using a bernoulli generated variable from numpy) then take the average of sequences that result at each flip, followed by using matplotlib to represent the findings in a graph? (Note the 'estimate of p' on the y-axis is an average of the sequence at each flip).

Upvotes: 0

Views: 251

Answers (1)

Nathan
Nathan

Reputation: 3648

You first create an array of 0s

x = np.zeros(nTrials+1, float)

and then fill that out.

x[i+1] =  np.random.binomial(1, p)

In the next sequence, your array is already filled out and you're overwritting data.

You have to

  1. create a new array for each sequence
  2. take the average over the simulated part of the array (leave out the trailing 0s)

You could also allow the mean list to not start at 0 to get a prettier graph.

Implementing

import numpy as np
import matplotlib.pyplot as plt
# Parameters

nSeq = 5
nTrials = 10 ** 3
p = 0.5  # Fair coin

# Containers

nrange = range(nTrials)

# Simulation

for j in range(nSeq):
    x = np.zeros(nTrials + 1, float)
    Mean_list = []  # re-initialize mean list for each sequence
    for i in range(nTrials):
        x[i] = np.random.binomial(1, p)
        xbar = np.mean(x[:i+1])
        Mean_list.append(xbar)

    plt.plot(nrange, Mean_list, label='Seq_' + str(j + 1))
plt.ylabel('Estimate of p')
plt.xlabel('Trials')
plt.legend(loc=0, ncol=3, fontsize='small')
plt.show()

Results

Upvotes: 2

Related Questions