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.