I am looking to find Pi using Monte Carlo Simulation. But somehow I got lost with my own code

Question

I am doing some online research and tutorial, then got stuck when i tried to find an approximation of pi using monte Carlo.

import math
import random
import time
import numpy as np
import matplotlib.pyplot as plt

random.seed(1234) 
old_est = 0
n = 1
MC_N= 1000000
results = np.repeat(0, MC_N)

for i in range(MC_N):
    x = random.random() 
    y = random.random()
    A_n = math.sqrt(x*2 + y*2)
    new_est = ((n-1)/n) * old_est + (1/n)*A_n
    results[n] = new_est
    if A_n <= 1:
        n +=1
    if n > MC_N:
        break
    old_est = new_est
    A_n = 4 * n / MC_N
print(results)
plt.plot(results)
plt.ylabel('numbers')
plt.show()

So I am using monte Carlo method and want to have an early stopping method to get the best estimation. I suppose to get around 3.14 but I got 0 Can somebody help me why I am wrong?

PS: It supposed to be calculating the approximation number of pi, To estimate the value of PI, we need the area of the square and the area of the circle. To find these areas, we will randomly place dots on the surface and count the dots that fall inside the circle and dots that fall inside the square. Such will give us an estimated amount of their areas. Therefore instead of using the actual areas, we will use the count of dots to use as areas. Also, I want to show that the error in estimation also decreased exponentially as the number of iterations increased.

how is the algorithm supposed to work? You haven't labeled anything — anon01, Apr 01 '21 at 03:58
pi, or phi?? Are you just counting points in circle vs square? What do your variables represent? — anon01, Apr 01 '21 at 04:06
You should step through your code and see what is happening. Basic debugging is your responsibility, not ours. — Mad Physicist, Apr 01 '21 at 06:01

score 0 · Answer 1 · answered Apr 01 '21 at 04:22

Here is a way to do it:

import numpy as np
import matplotlib.pyplot as plt
N = 1000000

# N random x, y points drawn uniformly from square(0,1)
x = np.random.rand(N)
y = np.random.rand(N)

# 1 if inside circle radius 1, 0 otherwise (first quadrant only)
circle_counts = (x**2 + y**2 < 1).astype(int)

# cumulative fraction of points that fall within circle region
pi_estimate = 4*circle_counts.cumsum()/np.arange(1,N+1)

# pi_estimate = array([4, 4, 4, ..., 3.14215828, 3.14215914, 3.14216])

plt.plot(x=np.arange(1,N+1), y=pi_estimate)
plt.show()

score 0 · Answer 2 · answered Apr 01 '21 at 04:55

I looked it up in wiki: https://en.wikipedia.org/wiki/Monte_Carlo_method It seems that your algorithm is very different, especially the line new_est = ((n-1)/n) * old_est + (1/n)*A_n. I am not sure whether this is what you want:

import math
import random
import time
import numpy as np
import matplotlib.pyplot as plt

random.seed(1234) 
MC_N = 100
results = []

est = 0
for i in range(1, MC_N):
    x = random.random() 
    y = random.random()
    A_n = math.sqrt(x**2 + y**2)    # not x*2 + y*2  
    if A_n < 1:
        est += 1    
        
    results.append(4 * est / MC_N)
    
print(results)
plt.plot(results)
plt.ylabel('numbers')
plt.show()

I am looking to find Pi using Monte Carlo Simulation. But somehow I got lost with my own code

2 Answers2