Strange increasing of sigma’s pi distribution increasing the time of repetition Python

I’m trying to perform this exercise 3.1 using python:

enter image description here

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:

enter image description here

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?

Answer

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:

enter image description here