Search code examples
pythonmatplotlibnumerical-methods

Using Monte Carlo method in python

I am working on the first version of the question written in the image below. I have generated a single random point using the rand command and tested whether or not that point was within the circle. Does my code accept the number of Monte Carlo functions as input value? I believe I have chosen N, the number of points to be small enough so I don't run out of memory. Also, when I run this code, it runs without any error but the graph does not show up. Looking for some help in places I might have went wrong.

enter image description here

import numpy as np
import matplotlib.pyplot as plt
from random import random

xinside = []
yinside = []
xoutside = []
youtside = []

insidecircle = 0
totalpoints = 10**3

for _ in range(totalpoints):
    x = random()
    y = random()
    if x**2+y**2 <= 1:
        insidecircle += 1
        xinside.append(x)
        yinside.append(y)
    else:
        xoutside.append(x)
        youtside.append(y)

fig, ax = plt.subplots()
ax.set_aspect('equal')
ax.scatter(xinside, yinside, color='g', marker='s')
ax.scatter(xoutside, youtside, color='r', marker='s')
fig.show()
Solution

  • The graph not showing is mysterious, perhaps try plt.show(). Alternatively, you can save the plot using savefig. Here is a working function (just modifying your posted code in the question) for the first part of your code together with the desired output plot.

    import numpy as np
import matplotlib.pyplot as plt
from random import random

def monte_carlo(n_points):
    xin, yin, xout, yout = [[] for _ in range(4)] # Defining all 4 lists together
    insidecircle = 0

    for _ in range(n_points):
        x = random()
        y = random()
        if x**2+y**2 <= 1:
            insidecircle += 1
            xin.append(x)
            yin.append(y)
        else:
            xout.append(x)
            yout.append(y)

    print ("The estimated value of Pi for N = %d is %.4f" %(n_points, 4*insidecircle/n_points))

    fig, ax = plt.subplots()
    ax.set_aspect('equal')
    ax.scatter(xin, yin, color='g', marker='o', s=4)
    ax.scatter(xout, yout, color='r', marker='o', s=4)
    plt.savefig('monte_carlo.png')

n_points = 10**4
monte_carlo(n_points)

> The estimated value of Pi for N = 10000 is 3.1380

    enter image description here

    Vectorized approach You can call it one-liner if you exclude the print statement in the function. I leave the time analysis as your homework

    import numpy as np
import matplotlib.pyplot as plt

def monte_carlo(n_points, x, y):
    pi_vec = 4*(x**2 + y**2 <= 1).sum()/n_points
    print ("The estimated value of Pi for N = %d is %.4f" %(n_points, pi_vec))

# Generate points
n_points = 10**4
x = np.random.random(n_points)
y = np.random.random(n_points)
# x = [random() for _ in range(n_points)] # alternative way to define x
# y = [random() for _ in range(n_points)] # alternative way to define y

# Call function
monte_carlo(n_points, x, y)

> The estimated value of Pi for N = 10000 is 3.1284

    Alternatively, you can get rid of two variables x and y by just using a single array of x and y points as follows:

    pts = np.random.uniform(0, 1, 2*n_points).reshape((n_points, 2))
monte_carlo(n_points, pts) # Function will be def monte_carlo(n_points, pts):

    And in the function use 

    pi_vec = 4*(pts[:,0]**2 + pts[:,1]**2 <= 1).sum()/n_points