CODEWITHSUNDEEP
pythonmontecarlo
Sookie
Sookie

Reputation: 29

Beginner Python Monte Carlo Simulation

I'm a beginner at Python and am working through exercises set by our instructor. I am struggling with this question.

In the Python editor, write a Monte Carlo simulation to estimate the value of the number π. Specifically, follow these steps: A. Produce two arrays, one called x, one called y, which contain 100 elements each, which are randomly and uniformly distributed real numbers between -1 and 1. B. Plot y versus x as dots in a plot. Label your axes accordingly. C. Write down a mathematical expression that defines which (x, y) pairs of data points are located in a circle with radius 1, centred on the (0, 0) origin of the graph. D. Use Boolean masks to identify the points inside the circle, and overplot them in a different colour and marker size on top of the data points you already plotted in B.

This is what I have at the moment.

import numpy as np
import math
import matplotlib.pyplot as plt
np.random.seed(12345)
x = np.random.uniform(-1,1,100) 
y = np.random.uniform(-1,1,100) 
plt.plot(x,y) //this works


for i in x:
    newarray = (1>math.sqrt(y[i]*y[i] + x[i]*x[i]))
plt.plot(newarray)

Any suggestions?

Upvotes: 3

Views: 2232

Answers (2)

jlandercy
jlandercy

Reputation: 11097

You are close to the solution. I slightly reshape your MCVE:

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

np.random.seed(12345)

N = 10000
x = np.random.uniform(-1, 1, N) 
y = np.random.uniform(-1, 1, N)

Now, we compute a criterion that makes sense in this context, such as the distance of points to the origin:

d = x**2 + y**2

Then we use Boolean Indexing to discriminate between points within and outside the Unit Circle:

q = (d <= 1)

At this point lies the Monte Carlo Hypothesis. We assume the ratio of uniformly distributed points in the Circle and in the plane U(-1,1)xU(-1,1) is representative for the Area of the Unit Circle and the Square. Then we can statistically assess pi = 4*(Ac/As) from the ratio of points within the Circle/Square. This leads to:

pi = 4*q.sum()/q.size # 3.1464

Finally we plot the result:

fig, axe = plt.subplots()
axe.plot(x[q], y[q], '.', color='green', label=r'$d \leq 1$')
axe.plot(x[~q], y[~q], '.', color='red', label=r'$d > 1$')
axe.set_aspect('equal')
axe.set_title(r'Monte Carlo: $\pi$ Estimation')
axe.set_xlabel('$x$')
axe.set_ylabel('$y$')
axe.legend(bbox_to_anchor=(1, 1), loc='upper left')
fig.savefig('MonteCarlo.png', dpi=120)

It outputs:

enter image description here

Upvotes: 4

Alessandro Solbiati
Alessandro Solbiati

Reputation: 979

as pointed out in the comment the error in your code is for i in x should be for i in xrange(len(x))

If you want to actually use a Boolean mask as said in the statement you could do something like this

    import pandas as pd
    allpoints = pd.DataFrame({'x':x, 'y':y})

    # this is your boolean mask
    mask = pow(allpoints.x, 2) + pow(allpoints.y, 2) < 1
    circlepoints = allpoints[mask]

    plt.scatter(allpoints.x, allpoints.y)
    plt.scatter(circlepoints.x, circlepoints.y)

increasing the number of point to 10000 you would get something like this scatter plot

to estimate PI you can use the famous montecarlo derivation

    >>> n = 10000
    >>> ( len(circlepoints) * 4 ) / float(n)
    <<< 3.1464

Upvotes: 6

Related Questions