This is the solution to Homework 10: Problems - Python advanced Monte Carlo.

The following figure illustrates the grade distribution for this homework.

Maximum possible points is 100.

This homework further explores Monte Carlo methods in Python.

1.  Monte Carlo approximation of the number $\pi$. Suppose we did not know the value of $\pi$ and we wanted to estimate its value using Monte Carlo methods. One practical approach is to draw a square of unit side, with its diagonal opposite corners extending from the coordinates origin $(0,0)$ to $(1,1)$. Now we try to simulate uniform random points from inside of this square by generating uniform random points along the $X$ and $Y$ axes, i.e., by generating two random uniform numbers (x,y) from the range $[0,1]$.

Now the generated random point $P$ has the coordinate $(x,y)$, so we can calculate its distance from the coordinate origin. Now suppose we also draw a quarter-circle inside of this square whose radius is unit and is centered at the origin $(0,0)$. The ratio of the area of this quarter-circle, $S_C$ to the area of the area of the square enclosing it, $S_S$ is,

This is because the area of the square of unit sides, is just 1. Therefore, if we can somehow measure the area of the quarter $S_C$, then we can use the following equation, to get an estimate of $\pi$,

In order to obtain, $S_C$, we are going to throw random points in the square, just as described above, and then find the fraction of points, $f=n_C/n_{\rm total}$, that fall inside this quarter-circle. This fraction is related to the area of the circle and square by the following equation,

Therefore, one can obtain an estimate of $\pi$ using this fraction,

Now, write a Python function, that takes in the number of points to be simulated, and then calculates an approximate value for $\pi$ based on the Monte Carlo algorithm described above. Write a second function that plot the estimate of $\pi$ versus the number of points simulated, like the following,

Here is an example Python code that estimates $\pi$.

import numpy.random as rnd
import numpy.linalg as nlg
import numpy as np
import matplotlib.pyplot as plt

def estimatePi(n=100000):
    x = np.zeros((n,4),dtype=np.dtype)
    counter = 0
    for i in range(n):
        x[i,0:2] = rnd.uniform(0.0,1.0,2)  # generate two uniform random numbers
        x[i,2] = nlg.norm(x[i,:], ord=2)
        if (x[i,2]<=1.0): counter += 1
        x[i,3] = 4.0*np.double(counter)/np.double(i+1)
    return x # the running approximate of pi is returned as a vector in the fourth column of x

def plotPi(n=10000):
    x = estimatePi(n)
    plt.semilogx( list(range(1,n+1)) \
                , x[:,3] \
                ) # plot with color red, as line
    plt.xlabel('Number of simulated points')
    plt.ylabel('Approximate value of Pi')
    #plt.axis([1, n , 0.1, 4.0]) # [xmin, xmax, ymin, ymax]
    plt.title('Estimating Pi by Monte Carlo simulation')


2.  Monte Carlo sampling of distribution functions Suppose that you wanted to generate points whose distribution follows the blue curve in the following curve, whose mathematical formulation is known (in the case here, the function is just the sum of two Gaussian distributions).

Now, one oway of doing this, is to draw a box around this curve, such that the box encompasses the entire curve.

Then, just as we did in the previous problem above, we draw random points from this square, and keep only those points that fall beneath this blue curve, like the red points in the following animation.

Now, if you plot the histogram of these points, you will see that the distribution of the red points follows closely the blue curve that we wanted to sample.

Now, given the above example, consider the following distribution function which we want to sample,

Suppose we know already that the highest point (maximum value) of this function is $f<0.2$, so that the value of this function always remains below $0.2$ everywhere along the positive x-axis, as seen in the following figure,

(A) Now, first write a function that generates a plot of this function, similar to the above plot.

Here is an example Python script that plots the requested curve.

import numpy as np
import matplotlib.pyplot as plt
x = np.linspace(0.001,15,100)
f= lambda x: np.exp(-(x-1)**2/2./x)*(x+1)/12.
fx = f(x)

(B) Then write another Python script, that samples from this function by first drawing a rectangle of base size $[0,15]$ and height $[0,h]$ with $h=0.2$. Then, draw uniform random points from this rectangle, and keep those that fall beneath the the value of $f(x)$ given above as points that are sampled from this function. Finally make a histogram of these points like the following.

Here is an example Python script that makes the requested plot.

u1 = np.random.rand(10000)*15 # uniform random samples scaled out
u2 = np.random.rand(10000)    # uniform random samples
idx = np.where(u2<=f(u1)/h)[0] # rejection criterion
v = u1[idx]

(C) Now make a plot of all generated points, both those that were accepted as samples, and those that were rejected, similar to the following plot, with accepted points in red color, and rejected points in black,

Here is an example Python script that makes the requested plot.