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 ther coordinate $(x,y)$, so we can calculate its distance from the coordiante 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,
\[\frac{S_C}{S_S} = \frac{\frac{1}{4}\pi r^2}{r^2} = \frac{1}{4}\pi\]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$,
\[\pi = 4S_C\]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 fracton is related to the area of the circle and square by the following equation,
\[f=\frac{n_C}{n_{\rm total}} = \frac{S_C}{S_S}\]Therefore, one can obtain an estimate of $\pi$ using this fraction,
\[\pi \approx \frac{1}{4}\frac{n_C}{n_{\rm total}}\]Now, write a Python function, that takes in the number of points to be simulated, and the calculate 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,
2. Monte Carlo sampling of distribution functions Suppose that you wanted to generate points whose distribution follows the blue curve in the following figure, 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,
\[f(x) = \frac{(x+1)}{12} \exp\bigg(-\frac{(x-1)^2}{2x}\bigg) ~~,~~ x > 0.\]Suppose we know already that the highest point (maximum value) of this function is $f_{\rm max}<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.
(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.
(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,
