In this notebook we calculate the value of $ \pi $ with a self coded random number generator.
In [28]:
import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline
%precision 8
plt.style.use('ggplot')
In [29]:
def rng(m=2**32, a=1103515245, c=12345):
rng.current = (a*rng.current + c) % m
return rng.current/m
# setting the seed
rng.current = 1
In [31]:
np.mean([np.exp(rng()) for i in range(1000000)])
Out[31]:
In [32]:
np.exp(1)
Out[32]:
In [33]:
np.var([np.exp(rng()) for i in range(1000000)])
Out[33]:
In [34]:
(np.exp(2)-1)*0.5-(np.exp(1)-1)**2
Out[34]:
In [35]:
import numpy as numpy
import matplotlib.pyplot as plt
import time as time
# input total number of random points and number of samples for the histogram
total_random_points = int(input("\nNumber of random points for Monte Carlo estimate of Pi?\n>"))
W = int(input("\nNumber of Samples for Monte Carlo estimate of Pi?\n>"))
# start time of calculation
start_time = time.time()
# number of random points inside unit cicrle and total random points
inside_circle = 0
# create empty x and y arrays for eventual scatter plot of generated random points
x_plot_array = numpy.empty(shape=(1,total_random_points))
y_plot_array = numpy.empty(shape=(1,total_random_points))
pi_approx = numpy.empty(shape=(1,total_random_points))
# generate random points and count points inside unit circle
# top right quadrant of unit cicrle only
for i in range(0, total_random_points):
# print(f'\nIteration: {i + 1}')
# generate random x, y in range [0, 1]
x = rng()
x_plot_array[0,i] = x
#numpy.append(x_plot_array, [x])
#print(f'{x_plot_array}')
# print(f'x value: {x}')
y = rng()
y_plot_array[0,i] = y
#numpy.append(y_plot_array, [y])
#print(f'{y_plot_array}')
# print(f'y value: {y}')
# calc x^2 and y^2 values
x_squared = x**2
y_squared = y**2
# count if inside unit circle, top right quadrant only
if numpy.sqrt(x_squared + y_squared) < 1.0:
inside_circle += 1
# print(f'Points inside circle {inside_circle}')
# print(f'Number of random points {i+1}')
# calc approximate pi value
pi_approx[0,i] = inside_circle / (i+1) * 4
#numpy.append(pi_approx,[inside_circle / (i+1) * 4])
# print(f'Approximate value for pi: {pi_approx}')
pi_approx_sample = numpy.empty(shape=(1,W))
for j in range(W):
x = numpy.array([rng() for i in range(total_random_points)])
y = numpy.array([rng() for i in range(total_random_points)])
x_squared = x**2
y_squared = y**2
pi_approx_sample[0,j] = 4*numpy.mean((x_squared+y_squared < 1.0))
# final numeric output for pi estimate
print ("\n--------------")
print (f"\nApproximate value for pi: {pi_approx[0,-1]}")
print (f"Difference to exact value of pi: {pi_approx[0,-1]-numpy.pi}")
print (f"Percent Error: (approx-exact)/exact*100: {(pi_approx[0,-1]-numpy.pi)/numpy.pi*100}%")
print (f"Execution Time: {time.time() - start_time} seconds\n")
# plot output of random points and circle, top right quadrant only
random_points_plot = plt.scatter(x_plot_array, y_plot_array, color='blue', s=.1)
circle_plot = plt.Circle( ( 0, 0 ), 1, color='red', linewidth=2, fill=False)
ax1 = plt.gca()
ax1.cla()
ax1.add_artist(random_points_plot)
ax1.add_artist(circle_plot)
plt.show()
plt.plot(pi_approx[0,1:]-numpy.pi)
plt.show()
plt.hist(pi_approx_sample[0,:]-numpy.pi,bins=int(numpy.sqrt(total_random_points)))
plt.show()
