Search code examples
pythonmontecarlo

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?

Solution

  • 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