Imagine that we want to measure the area of a pond with arbitrary shape. Suppose that this pond is in the middle of a field with known area $A$. If we throw $N$ stones randomly, such that they land within the boundaries of the field, and we count the number of stones that fall in the pond $N_{in}$, the area of the pond will be approximately proportional to the fraction of stones that make a splash, multiplied by $A$: $$A_{pond}=\frac{N_{in}}{N}A.$$ This simple procedure is an example of the “Monte Carlo” method.
More generaly, imagine a rectangle of height $H$ in the integration interval $[a,b]$, such that the function $f(x)$ is within its boundaries. Compute $n$ pairs of random numbers $(x_i,y_i)$ such that they are uniformly distributed inside this rectangle. The fraction of points that fall within the area contained below $f(x)$, *i. e.*, that satisfy $y_i \leq f(x_i)$ is an estimate of the ratio o fthe integral of $f(x)$ and the area of the rectangle. Hence, the estimate of the integral will be given by: $$\int _a^b{f(x)dx} \simeq I(N) = \frac{N_{in}}{N}H(b-a). $$
Another Monte Carlo procedure is based on the definition: $$\langle g \rangle=\frac{1}{(b-a)} \int _a^b{f(x)dx}. $$ In order to determine this average, we sample the value of $f(x)$: $$\langle f \rangle \simeq \frac{1}{N}\sum_{i=1}^{N}f(x_i),$$ where the $N$ values $x_i$ are distributed unformly in the interval $[a,b]$. The integral will be given by $$I(N)=(b-a) \langle f \rangle .$$
The Monte Carlo method clearly yields approximate results. The accuracy deppends on the number of values $N$ that we use for the average. A possible measure of the error is the “variance” $\sigma^2$ defined by: $$\sigma ^2=\langle f^2 \rangle - \langle f \rangle ^2, $$ where $$\langle f \rangle = \frac{1}{N} \sum_{i=1}^N f(x_i)$$ and $$\langle f^2 \rangle = \frac{1}{N} \sum_{i=1}^{N} f(x_i)^2.$$ The “standard deviation” is $\sigma$. However, we should expect that the error decreases with the number of points $N$, and the quantity $\sigma$ defines by ([mc_sigma]) does not. Hence, this cannot be a good measure of the error.
Imagine that we perform several measurements of the integral, each of them yielding a result $I_n$. These values have been obtained with different sequences of $N$ random numbers. According to the central limit theorem, these values whould be normally dstributed around a mean $\langle I \rangle$. Suppouse that we have a set of $M$ of such measurements ${I_n}$. A convenient measure of the differences of these measurements is the “standard deviation of the means” $\sigma_M$: $$\sigma_M ^2=\langle I^2 \rangle - \langle I \rangle ^2, $$ where $$\langle I \rangle = \frac{1}{M} \sum_{n=1}^M I_n$$ and $$\langle I^2 \rangle = \frac{1}{M} \sum_{n=1}^{M} I_n^2.$$ It can be proven that $$\sigma_M \approx \sigma/\sqrt{N}. $$ This relation becomes exact in the limit of a very large number of measurements. Note that this expression implies that the error decreases with the square root of the number of trials, meaning that if we want to reduce the error by a factor 10, we need 100 times more points for the average.
Write a program that implements the “hit and miss” Monte Carlo integration algorithm. Find the estimate $I(N)$ for the integral of $$f(x)=4\sqrt{1-x^2}$$ as a function of $N$, in the interval $(0,1)$. Choose $H=1$, and sample only the $x$-dependent part $\sqrt{1-x^2}$, and multiply the result by 4. Calculate the difference between $I(N)$ and the exact result $\pi$. This difference is a measure of the error associated with the Monte Carlo estimate. Make a log-log plot of the error as a function of $N$. What is the approximate functional deppendece of the error on $N$ for large $N$?
Estimate the integral of $f(x)$ using the simple Monte Carlo integration by averaging over $N$ points, using ([mc_integral2]), and compute the error as a function of $N$, for $N$ upt to 10,000. Determine the approximate functional deppendence of the error on $N$ for large $N$. How many trials are necessary to determine $I_N$ to two decimal places?
Perform 10 measurements $I_n(N)$, with $N=10,000$ using different random sequences. Show in a table the values of $I_n$ and $\sigma$ according to ([mc_integral2]) and ([mc_sigma]). Use ([mc_sigmam]) to estimate the standard deviation of the means, and compare to the values obtained from ([mc_sigma2]) using the 100,000 values.
To verify that your result for the error is independent of the number of sets you used to divide your data, repeat the previous item grouping your results in 20 groups of 5,000 points each.
To examine the effects of a poor random number generator, modify your program to use the linear congruential random number generator using the perameters $a=5$, $c=0$ and the seed $x_1=1$. Repeat the integral of the previous exercise and compare your results.
Exercise 10.2
In [1]:
%matplotlib inline
import numpy as np
from matplotlib import pyplot
x = np.arange(0,1,0.02)
pyplot.plot(x, 4*np.sqrt(1-x**2))
Out[1]:
In [2]:
# Hit and miss Monte Carlo integration
ngroups = 16
I = np.zeros(ngroups)
N = np.zeros(ngroups)
E = np.zeros(ngroups)
n0 = 100
for i in range(ngroups):
N[i] = n0
x = np.random.random(n0)
y = np.random.random(n0)
I[i] = 0.
Nin = 0
for j in range(n0):
if(y[j] < np.sqrt(1-x[j]**2)):
Nin += 1
I[i] = 4.*float(Nin)/float(n0)
E[i] = abs(I[i]-np.pi)
print (n0,Nin,I[i],E[i])
n0 *= 2
pyplot.plot(N,E,ls='-',c='red',lw=3);
pyplot.plot(N,0.8/np.sqrt(N),ls='-',c='blue',lw=3);
pyplot.xscale('log')
pyplot.yscale('log')
In [3]:
# Simple Monte Carlo Integration
ngroups = 16
I = np.zeros(ngroups)
N = np.zeros(ngroups)
E = np.zeros(ngroups)
n0 = 100
for i in range(ngroups):
N[i] = n0
r = np.random.random(n0)
I[i] = 0.
for j in range(n0):
x = r[j]
I[i] += np.sqrt(1-x**2)
I[i] *= 4./float(n0)
E[i] = abs(I[i]-np.pi)
print (n0,I[i],E[i])
n0 *= 2
pyplot.plot(N,E,ls='-',c='red',lw=3);
pyplot.plot(N,0.8/np.sqrt(N),ls='-',c='blue',lw=3);
pyplot.xscale('log')
pyplot.yscale('log')
In [5]:
n0 = 100000
I = np.zeros(n0)
r = np.random.random(n0)
for j in range(n0):
x = r[j]
I[j] = 4*np.sqrt(1-x**2)
def group_measurements(ngroups):
global I,n0
nmeasurements = n0//ngroups
for n in range(ngroups):
Ig = 0.
Ig2 = 0.
for i in range(n*nmeasurements,(n+1)*nmeasurements):
Ig += I[i]
Ig2 += I[i]**2
Ig /= nmeasurements
Ig2 /= nmeasurements
sigma = Ig2-Ig**2
print(Ig,Ig2,sigma)
group_measurements(10)
print("=============================")
group_measurements(20)
print("=============================")
group_measurements(1)
If the function being integrated does not fluctuate too much in the interval of integration, and does not differ much from the average value, then the standard Monte Carlo mean-value method should work well with a reasonable number of points. Otherwise, we will find that the variance is very large, meaning that some points will make small contributions, while others will make large contributions to the integral. If this is the case, the algorithm will be very inefficient. The method can be improved by splitting the function $f(x)$ in two $f(x)=f_1(x)+f_2(x)$, such that the integral of $f_1(x)$ is known, and $f_2(x)$ as a small variance. The “variance reduction” technique, consists then in evaluating the integral of $f_2(x)$ to obtain: $$\int _a^b{f(x)dx}=\int _a^b {f_1(x)dx} + \int _a^b{f_2(x)dx} = \int _a^b{f_1(x)dx}+J.$$
Imagine that we want to sample the function $f(x)=e^{-x^2}$ in the interval $[0,1]$. It is evident that most of our points will fall in the region where the value of $f(x)$ is very small, and therefore we will need a large number of values to achieve a decent accuracy. A way to improve the measurement by reducing the variance is obtained by “importance sampling”. As the name says, the idea is to sample the regions with larger contributions to the integral. For this goal, we introduce a probability distribution $P(x)$ normalized in the interval of integration $$\int _a^b{P(x)dx} = 1.$$ Then, we can rewrite the integral of $f(x)$ as $$I=\int _a^b{\frac{f(x)}{P(x)}P(x)dx} $$ We can evaluate this integral, by sampling according to the probability distribution $P(x)$ and evaluating the sum $$I(N)=\frac{1}{N} \sum_{i=1}^N \frac{f(x_i)}{P(x_i)}. $$ Note that for the uniform case $P(x)=1/(b-a)$, the expression reduces to the simple Monte Carlo integral.
We are free to choose $P(x)$ now. We wish to do it in a way to reduce and minimize the variance of the integrand $f(x)/P(x)$. The way to to this is picking a $P(x)$ that mimics $f(x)$ where $f(x)$ is large. if we are able to determine an apropiate $P(x)$, the integrand will be slowly varying, and hence the variance will be reduced. Another consideration is that the generation of points according to the distribution $P(x)$ should be a simple task. As an example, let us consider again the integral $$I=\int _0^1 {e^{-x^2}dx}.$$ A reasonable choice for a weigh function is $P(x)=Ae^{-x}$, where $A$ is a normalization constant.
Notice that for $P(x)=f(x)$ the variance is zero! This is known as the zero variance property. There is a catch, though: The probability function $P(x)$ needs to be normalized, implying that in reality, $P(x)=f(x)/\int f(x)dx$, which assumes that we know in advance precisely the integral that we are trying to calculate!
Choose the weight function $P(x)=e^{-x}$ and evaluate the integral: $$\int _0^{\infty} {x^{3/2}e^{-x}dx}.$$
Choose $P(x)=e^{-ax}$ and estimate the integral $$\int _0^{\pi} \frac{dx}{x^2+\cos ^2{x}}.$$ Determine the value of $a$ that minimizes the variance of the integral.
In [16]:
pyplot.xlim(0,10)
x = np.arange(0,10,0.1)
pyplot.plot(x,np.exp(-x));
pyplot.plot(x,np.exp(-x**2));
pyplot.plot(x,x**1.5*np.exp(-x));
In [9]:
# Trapezoidal integration
def trapezoids(func, xmin, xmax, nmax):
Isim = func(xmin)+func(xmax)
h = (xmax-xmin)/nmax
for i in range(1,nmax):
x = xmin+i*h
Isim += 2*func(x)
Isim *= h/2
return Isim
def f(x):
return x**1.5*np.exp(-x)
print("Trapezoids: ", trapezoids(f, 0., 20., 100000))
# Simple Monte Carlo integration
n0 = 1000000
r = np.random.random(n0)
Itot = np.sum(r**1.5*np.exp(-r))
print("Simple Monte Carlo: ", Itot/n0)
x = -np.log(r)
Itot = np.sum(x**1.5)
print("Importance Sampling: ", Itot/n0)
In [11]:
pyplot.xlim(0,np.pi)
x = np.arange(0,np.pi,0.05)
pyplot.plot(x,1./(x**2+np.cos(x)**2));
pyplot.plot(x,np.exp(-x));
pyplot.plot(x,np.exp(-2*x));
pyplot.plot(x,np.exp(-0.2*x));
In [8]:
# Trapezoidal integration
def g(x):
return 1./(x**2+np.cos(x)**2)
print("Trapezoids: ", trapezoids(g, 0., np.pi, 1000000))
# Simple Monte Carlo integration
n0 = 1000000
a = np.arange(0.1,2.1,0.1)
I = np.arange(0.1,2.1,0.1)
r = np.random.random(n0)
I0 = np.sum(1./((r*np.pi)**2+np.cos(r*np.pi)**2))
print("Simple Monte Carlo: ", I0/n0*np.pi)
# Importance Sampling
print("Importance Sampling:")
x = -np.log(r)
i = 0
for ai in a:
norm = (1.-np.exp(-ai*np.pi))/ai
x1 = norm*x/ai
Itot = 0.
Nin = 0
I2 = 0.
for xi in x1:
if(xi <= np.pi):
Nin += 1
Itot += g(xi)*np.exp(xi*ai)
I2 += (g(xi)*np.exp(xi*ai))**2
Itot *= norm
I2 *= norm
I[i] = Itot/Nin
i += 1
print(ai,Itot/Nin,np.sqrt(abs(Itot**2/Nin**2-I2/Nin))/np.sqrt(Nin))
pyplot.plot(a,I,ls='-',marker='o',c='red',lw=3);
Use the Metropolis algorithm to sample points according to a ditribution and estimate the integral $$\int _0^4 {x^2e^{-x}dx},$$ with $P(x)=e^{-x}$ for $0 \leq x \leq 4$. Plot the number of times the walker is at points $x_0$, $x_1$, $x_2$, ... Is the integrand sampled uniformly? If not, what is the approximate region of $x$ where the integrand is sampled more often?
In [26]:
delta = 2
xmin = 0.
xmax = 4.
def f(x):
return x**2*np.exp(-x)
def P(x):
global xmin, xmax
if(x < xmin or x > xmax):
return 0.
return np.exp(-x)
def metropolis(xold):
global delta
xtrial = np.random.random()
xtrial = xold+(2*xtrial-1)*delta
weight = P(xtrial)/P(xold)
xnew = xold
if(weight >= 1): #Accept
xnew = xtrial
elif(weight != 0):
r = np.random.random()
if(r <= weight): #Accept
xnew = xtrial
return xnew
xwalker = (xmax+xmin)/2.
for i in range(100000):
xwalker = metropolis(xwalker)
I0 = 0.
N = 300000
x = np.zeros(N)
x[0] = xwalker
for i in range(1,N):
for j in range(20):
xwalker = metropolis(xwalker)
x[i] = xwalker
I0 += x[i]**2
binwidth=0.1
pyplot.hist(x,bins=np.arange(xmin-1, xmax+1, 0.1),normed=True);
print("Trapezoids: ", trapezoids(f,xmin,xmax,100000))
print("Metropolis: ", I0*(1.-np.exp(-4.))/N)
In [23]:
pyplot.hist(x**2,bins=np.arange(xmin**2-1, xmax**2+1, 0.1),normed=True);
In [ ]: