pythongaussianpi
Strange increasing of sigma's pi distribution increasing the time of repetition Python
I'm trying to perform this exercise 3.1 using python:
the code "works" and is the following:
#NUMERICAL ESTIMATE OF PI
import numpy as np #library for numerical calculations
import matplotlib.pyplot as plt #library for plotting purposes
from scipy.stats import norm #needed for gaussian fit
#*******************************************************************************
M = 10**2 #number of times we calculate pi
N = 10**4 #number of point generated
mean_pi=[] #empy list
for i in range(M): #for loops over the period
x=np.random.uniform(-1,1,N) #array of the given shape d and populate it with random samples from a uniform distribution over [-1,1)
y=np.random.uniform(-1,1,N) #array of the given shape d and populate it with random samples from a uniform distribution over [-1,1)
x_sel=x[(x**2+y**2)<=1] #selection of x point
y_sel=y[(x**2+y**2)<=1] #selection of y point
mean_pi+=[4*len(x_sel)/len(x)] #list of pi's mean value
#*******************************************************************************
plt.figure(figsize=(8,3)) #a unique identifier for the figure
_,bins,_=plt.hist(mean_pi,bins=int(np.sqrt(N)),density=True, color="skyblue") #sintex to create a histogram from a dataset x with n bins
#and store an array specifying the bin ranges in the variable bins.
mu, sigma = norm.fit(mean_pi) #get the mean and standard deviation of data
k = sigma*np.sqrt(N) #k parameters
best_fit_line = norm.pdf(bins, mu, sigma) #get a line of best fit for the data
print("\nTime of repetitions:", M, ". The mean of the distribution is: ", mu, ". The standard deviation is:", sigma, ". The k parameters is:", k ,". \n")
#*******************************************************************************
plt.plot(bins, best_fit_line, color="red") #plot y versus x as lines and/or markers
plt.grid() #configure the grid lines
plt.xlabel('Bins',fontweight='bold') #set the label for the x-axis
plt.ylabel('Pi',fontweight='bold') #set the label for the y-axis
plt.title("Histogram for Pi vs. bins") #set a title for the scatter plot
plt.show() #display all open figures
print("\n")
#*******************************************************************************
M = 10**3 #number of times we calculate pi
N = 10**4 #number of point generated
mean_pi=[] #empy list
for i in range(M): #for loops over the period
x=np.random.uniform(-1,1,N) #array of the given shape d and populate it with random samples from a uniform distribution over [-1,1)
y=np.random.uniform(-1,1,N) #array of the given shape d and populate it with random samples from a uniform distribution over [-1,1)
x_sel=x[(x**2+y**2)<=1] #selection of x point
y_sel=y[(x**2+y**2)<=1] #selection of y point
mean_pi+=[4*len(x_sel)/len(x)] #list of pi's mean value
#*******************************************************************************
plt.figure(figsize=(8,3)) #a unique identifier for the figure
_,bins,_=plt.hist(mean_pi,bins=int(np.sqrt(N)),density=True, color="skyblue") #sintex to create a histogram from a dataset x with n bins
#and store an array specifying the bin ranges in the variable bins.
mu, sigma = norm.fit(mean_pi) #get the mean and standard deviation of data
k = sigma*np.sqrt(N) #k parameters
best_fit_line = norm.pdf(bins, mu, sigma) #get a line of best fit for the data
print("Time of repetitions:", M, ". The mean of the distribution is: ", mu, ". The standard deviation is:", sigma, ". The k parameters is:", k ,". \n")
#*******************************************************************************
plt.plot(bins, best_fit_line, color="red") #plot y versus x as lines and/or markers
plt.grid() #configure the grid lines
plt.xlabel('Bins',fontweight='bold') #set the label for the x-axis
plt.ylabel('Pi',fontweight='bold') #set the label for the y-axis
plt.title("Histogram for Pi vs. bins") #set a title for the scatter plot
plt.show() #display all open figures
print("\n")
#*******************************************************************************
M = 5*10**3 #number of times we calculate pi
N = 10**4 #number of point generated
mean_pi=[] #empy list
for i in range(M): #for loops over the period
x=np.random.uniform(-1,1,N) #array of the given shape d and populate it with random samples from a uniform distribution over [-1,1)
y=np.random.uniform(-1,1,N) #array of the given shape d and populate it with random samples from a uniform distribution over [-1,1)
x_sel=x[(x**2+y**2)<=1] #selection of x point
y_sel=y[(x**2+y**2)<=1] #selection of y point
mean_pi+=[4*len(x_sel)/len(x)] #list of pi's mean value
#*******************************************************************************
plt.figure(figsize=(8,3)) #a unique identifier for the figure
_,bins,_=plt.hist(mean_pi,bins=int(np.sqrt(N)),density=True, color="skyblue") #sintex to create a histogram from a dataset x with n bins
#and store an array specifying the bin ranges in the variable bins.
mu, sigma = norm.fit(mean_pi) #get the mean and standard deviation of data
k = sigma*np.sqrt(N) #k parameters
best_fit_line = norm.pdf(bins, mu, sigma) #get a line of best fit for the data
print("Time of repetitions:", M, ". The mean of the distribution is: ", mu, ". The standard deviation is:", sigma, ". The k parameters is:", k ,". \n")
#*******************************************************************************
plt.plot(bins, best_fit_line, color="red") #plot y versus x as lines and/or markers
plt.grid() #configure the grid lines
plt.xlabel('Bins',fontweight='bold') #set the label for the x-axis
plt.ylabel('Pi',fontweight='bold') #set the label for the y-axis
plt.title("Histogram for Pi vs. bins") #set a title for the scatter plot
plt.show() #display all open figures
#*******************************************************************************
print("\n How many couples N you need to estimate pi at better than 0.0001? The number of couples N is:", (k**2)*10**8 ,".")
#*******************************************************************************
With the output:
As you can see, the sigma increase, meanwhile i expect that it decrease when the time repetition increase...i don't understand where is the error.
I also tryed to increase N but the results are not better...
Someone can help me please?
Solution
I understand the error. In order to implement the right code it is necessary to chage the line:
_,bins,_=plt.hist(mean_pi,bins=int(np.sqrt(N)),density=True, color="slateblue") #sintex to create a histogram from a dataset x with n bins
in:
_,bins,_=plt.hist(mean_pi,bins=int(np.sqrt(M)),density=True, color="slateblue") #sintex to create a histogram from a dataset x with n bins
and the output will be:
