November 21, 2006---Class 23---Monte Carlo Error Analysis
Activities:
Monte Carlo Error Analysis (Section 11.4)
-----------------------------------------
In the past, I have assigned a homework problem related to this;
however, since I don't want to burden the class with more
homework, I will demonstrate how problem 11.10 on page 428 is
done.
(a) estimate the integral of exp(-x) between x=0 and 1 using the
sample mean Monte Carlo method. Use n=10^4, 10^6 and 10^8
points.
A code exp.c can be found on Nations in ~sg/nov21. This code
was explained in detail. The code was run with n=10^6 and 10^4
points, but 10^8 was taking too long. Since the program has an
optional second argument to repeat the experiment multiple
times, we generated 100 samples for each value of n. We also
run with n=10^3 since 10^8 was not available.
The value of the integral was computed using Mathematica.
It is 0.632121. The code computes both the estimate of the
integral and the variance of the function.
We plotted histograms of the results for all three values of n
using a simple shell script makehist.csh. We modified the
histogram to plot the error rather than the value of the
integral. Here are the commands to plot the histogram of error:
#!/bin/csh
(awk '{print $2-.632121}' exp_1000 | hist -n 20 -g 2; \
awk '{print $2-.632121}' exp_10000 | hist -n 20 -g 2 -m 2; \
awk '{print $2-.632121}' exp_1000000 | hist -n 20 -g 2 -m 3) |axis |plot -Tx
We found that the variance of the function divided by sqrt(n)
[See Eq 11.23b] was in very good agreement with the widths of
the three histograms. Thus we verified that we understand the
error from this method. It only falls like 1/sqrt(n), and the
error is smaller for functions with smaller variance. This
leads to the idea of important sampling methods discussed in
Section 11.6.
(e) We then looked at the integral of cos^2(theta) for theta
between 0 and 2 Pi. This integral is clearly 0.5. The book
asks us to do the integral by using a series of correlated
values for theta, i.e., a random walk in theta. When we do
this, we find the naive error formula Eq. 11.23 is not appropriate.
We looked at the codes cos_sq_ralk.c and cos_sq.c. The former
has the random walk, while the latter is a standard sample mean
Monte Carlo code. We ran cos_sq_rwalk with n=10^6 for 100
samples.
The code details.c prints the entire history of the random walk
and it also does at the same time a uniform sampling such as in
cos_sq.c. On each line of output is the value for each
procedure with the random walk first.
I did a long run, that filled my disk space. Unfortunately,
when the output was analyzed with blockerr, it was clear that
the error was not correctly calculated. Here is an example from
the random walk (first column) with 100000 samples:
100000 items, block size = 1 average= 4.948896e-01 +/- 0.000000e+00
This was not a good way to break for Thanksgiving. The problem
is due to round off error and will be fixed in the next class.