import numpy as np
import matplotlib.pyplot as plt
import random as rand
import statistics
import math
from scipy.optimize import curve_fit
1(a)
def sample(n):
return [rand.randrange (-1,2,2) for _ in range(n)]
#return n samples that are +1 or -1 with equal probability
#print([rand.randrange(-1,2,2) for _ in range(10)])
# for in range: repeat n times
#print([rand.randrange (-1,2,2) for _ in range(10)])
print (sample(10))
1(b)
# Define x
def meanx(n):
return sum(sample(n))/n
samplelist10=[]
# n=10
for i in range(1000):
samplelist10.append(meanx(10))
#n=100
samplelist100=[]
for i in range(1000):
samplelist100.append(meanx(100))
#n=1000
samplelist1000=[]
for i in range(1000):
samplelist1000.append(meanx(1000))
#n=10000
samplelist10000=[]
for i in range(1000):
samplelist10000.append(meanx(10000))
# Plot histograms
fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(nrows=2,ncols=2)
# This is the simplest way to create a plot with Pyplot. The 'fig' object represents the entire figure you will
# create, and the 'ax' object (a tuple of tuples here since we have two rows and two columns) represents the axes you will plot on.
fig.suptitle('Histograms for Various n')
# give specifics for each subplot
ax1.set_title('10 Samples')
ax1.set_xlim([-1, 1]) # Set the x-axis view limits
ax1.hist(samplelist10, density=True) # Make a histogram
# similarly for ax2, ax3 and ax4
ax2.set_title('100 Samples')
ax2.set_xlim([-1, 1]) # Set the x-axis view limits
ax2.hist(samplelist100, density=True) # Make a histogram
ax3.set_title('1000 Samples')
ax3.set_xlim([-1, 1]) # Set the x-axis view limits
ax3.hist(samplelist1000, density=True) # Make a histogram
ax4.set_title('10000 Samples')
ax4.set_xlim([-1, 1]) # Set the x-axis view limits
ax4.hist(samplelist10000, density=True) # Make a histogram
plt.tight_layout()
plt.show()
1(c)
#Compute varianace
var10=statistics.variance(samplelist10)
var100=statistics.variance(samplelist100)
var1000=statistics.variance(samplelist1000)
var10000=statistics.variance(samplelist10000)
# A minimal implementation of pyplot
#n=10
plt.plot(math.log(10), math.log(var10), 'o')
#n=100
plt.plot(math.log(100), math.log(var100), '*')
#n=1000
plt.plot(math.log(1000), math.log(var1000), 'D')
#n=10000
plt.plot(math.log(10000), math.log(var10000), '+')
plt.legend(['n=10','n=100','n=1000','n=10000'])
plt.xlabel(r'$log_{10} (n)$')
plt.ylabel(r'$log_{10} (var(\bar{x}))$')
plt.show()
1(d)
#Linear fit
# Plot the original data points
#n=10
plt.plot(math.log(10), math.log(var10), 'o')
#n=100
plt.plot(math.log(100), math.log(var100), '*')
#n=1000
plt.plot(math.log(1000), math.log(var1000), 'D')
#n=10000
plt.plot(math.log(10000), math.log(var10000), '+')
plt.legend(['n=10','n=100','n=1000','n=10000'])
plt.xlabel(r'$log_{10} (n)$')
plt.ylabel(r'$log_{10} (var(\bar{x}))$')
# Plot the fit function here
x=[math.log(10),math.log(100),math.log(1000),math.log(10000)]
y=[math.log(var10),math.log(var100),math.log(var1000),math.log(var10000)]
m, b = np.polyfit(x, y, 1)
def f(a):
return m*a+b
print(m)
i=0
pf=[]
for i in range(len(x)):
pf.append(f(x[i]))
plt.plot(x,pf)
plt.show()
From above, slope=alpha=m~-1 n^alpha~n^(-1)=1/n As n increases, n^alpha would decrease. Therefore, Var(xn) would decrease as n increasing.
1(e)
The unit of p(x) is 1/unit of x, since probability is unitless. When extract a probability from a histogram, the probability is determined as the area under the histogram. i.e, probability is the integral of p(x)dx over the the range x to x+dx. Therefore, if the bin width is small, we can estimate the probability by multiplying p(x) by the bin width.
1(f)
#n=10
h10,bin10,p10=ax1.hist(samplelist10,density=True) #h10: heights, bin10:list of edge
density10=[]
for i in range(len(h10)):
if h10[i]==0:
density10.append(0)
else:
density10.append(pow(10,m)*(math.log(h10[i])))
print(hy10)
#n=100
h100,bin100,p100=ax2.hist(samplelist100,density=True)
density100=[]
for i in range(len(h100)):
if h100[i]==0:
density100.append(0)
else:
density100.append((pow(100,m))*(math.log(h100[i])))
#n=1000
h1000,bin1000,p1000=ax3.hist(samplelist1000,density=True)
density1000=[]
for i in range(len(h1000)):
if h1000[i]==0:
density1000.append(0)
else:
density1000.append((pow(1000,m))*(math.log(h1000[i])))
#n=10000
h10000,bin10000,p10000=ax4.hist(samplelist10000,density=True)
density10000=[]
for i in range(len(h10000)):
if h10000[i]==0:
density10000.append(0)
else:
density10000.append((pow(10000,m))*(math.log(h10000[i])))
#-x^2/2
def ref(x):
return -(x*x)/2
refy=[]
for i in range(len(x)):
refy.append(ref(x[i]))
plt.xlabel(r"$\bar{x}_n$")
plt.ylabel(r"$n log(p(\bar{x}_n))$")
bin10=np.delete(bin10,-1)
plt.plot(bin10, density10, label='n=10')
bin100=np.delete(bin100,-1)
plt.plot(bin100, density100, label='n=100')
bin1000=np.delete(bin1000,-1)
plt.plot(bin1000, density1000, label='n=1000')
bin10000=np.delete(bin10000,-1)
plt.plot(bin10000, density10000,label='n=10000')
plt.plot(x,refy)
plt.legend(['n=10','n=100','n=1000','n=10000','$-x^2/2$'])
plt.show()
Question 2
2a:
def pd1(x):
return (1/(math.sqrt(2*(math.pi)*0.1)))*(math.exp(-(x*x)/(2*0.1)))
#n=10
x=np.linspace(-1,1,10000)
pd10=[]
for i in range(len(x)):
pd10.append(pd1(x[i]))
#n=100
def pd2(x):
return (1/(math.sqrt(2*(math.pi)*0.01)))*(math.exp(-x*x/(2*0.01)))
pd100=[]
for i in range(len(x)):
pd100.append(pd2(x[i]))
#n=1000
def pd3(x):
return (1/(math.sqrt(2*(math.pi)*0.001)))*(math.exp(-x*x/(2*0.001)))
pd1000=[]
for i in range(len(x)):
pd1000.append(pd3(x[i]))
#n=10000
def pd4(x):
return (1/(math.sqrt(2*(math.pi)*0.0001)))*(math.exp(-x*x/(2*0.0001)))
pd10000=[]
for i in range(len(x)):
pd10000.append(pd4(x[i]))
plt.xlabel("x")
plt.ylabel("Probability density function$")
plt.plot(x, pd10, label='n=10')
plt.plot(x,pd100, label='n=100')
plt.plot(x,pd1000, label='n=1000')
plt.plot(x,pd10000, label='n=10000')
plt.legend(['n=10','n=100','n=1000','n=10000'])
plt.show()
This plot is consistent with the histograms obtained in 1(b).The peak becomes sharper as n becomes larger, while the peak location stays constant.