MAT4376-HW1
# Code to calculate sample mean
print(np.mean(ntrials))
2.387
Appendix - Code
#Import libraries
import numpy as np
import numpy.random as npr
import matplotlib.pyplot as plt
from scipy.stats import beta
import matplotlib.pyplot as plt
import numpy as np
# plotting beta distribution (3,5)
a = 3
b = 5
x = np.linspace(beta.ppf(0.01, a, b), beta.ppf(0.99, a, b), 100)
plt.figure(figsize=(6, 3))
plt.title("Beta Distribution cdf - B(3,5)")
plt.plot(x, beta.pdf(x, a, b))
# rejection algorithm
n = 2000
M = 2.305
x = np.zeros(n)
ntrials = np.zeros(n)
for i in range(n):
flag = True # flag changes to False when distribution condition is satisfied
count = 0 # count the number of trials until condition is satisfied
while(flag):
u = npr.random()
y = npr.random()
if (u < (105*pow(y,2)*pow(1-y,4))/M) :
flag = False
count = count + 1
x[i] = y
ntrials[i] = count
# histogram for simulated sample
plt.figure(figsize=(6, 3))
plt.hist(x,bins=25,density=True,rwidth=.9)
s=np.arange(start=min(x),step=.01,stop=max(x))
## density of the standard Gaussian distribution
plt.plot(s,105*pow(s,2)*pow(1-s,4),'r', label='density of \n B(3,5)')
plt.legend()
plt.show()
# histogram for simulated sample trials
plt.figure(figsize=(6, 3))
plt.hist(ntrials,bins=np.arange(start=.5,step=1,stop=max(ntrials)+1),
density=True,rwidth=.9)
# theoretical pmf of the number of trials
M = 2.305
p=1/M
nn=np.arange(start=1,step=1,stop=max(ntrials)+1)
plt.plot(nn,p*(1-p)**(nn-1),'ro', label=r'pmf of the geometric distribution with mean $M$')
plt.legend()