CODEWITHSUNDEEP
pythonmatplotlibnumerical-methods
jh123
jh123

Reputation: 173

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()

Upvotes: 3

Views: 1318

Answers (1)

Sheldore
Sheldore

Reputation: 39042

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

Upvotes: 3

Related Questions