# Code to calculate sample mean print(np.mean(ntrials))

#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()