November 8, 2005 --- Class 20 --- Simpson's Rule Rederived, Hit or Miss Method

Activities:
	Simpson's Rule Rederived

	In the last class, starting from the trapezoid rule, we used 
	the formulae below for correcting the leading 1/N^2 errors.

	S_N  = A + E/N^2
	S_2N = A + E/(2N)^2

	so   A = 4/3 S_2N -1/3 S_N

	Applying this to a one-dimensional itegral, we easily found that 
	the weighting of the sum of two expressions on the right hand side
	exactly coincides with the Simpson's rule formula.  We previously
	found that Simpson's Rule has a leading error of order 1/N^4.  Since
	the approach above removes the leading 1/N^2 error and since the
	error of the trapezoid rule only has even powers of 1/N (asserted, but
	not proved), this is not surprising.

	Hit-or-Miss Method

	The Hit-or-Miss method and Sample Mean methods are described in 
	Section 11.2.  Let's try the first method and see how well it does.

	Program ~sg/src/pi1.c implements the hit-or-miss method using a
	random number generator that was in the first edition of Numerical
	Recipes.  It is called ran1 and the source is in ~sg/src/ran1.c.
	You have to tell the program how many "darts you want to throw."
	This is the first command line argument.  You may also have it
	repeat the experiment multiple times.  This is the optional second
	argument.  Look at the code to see how this is done if you are
	unfamiliar with this type of feature.  Also note how the random
	number generator gets its seed through a call to pid().  This is
	the "process id" of the calling routine.  It is essentially a
	random number.  If you like, you can modify the code so that you
	specify the seed yourself.  This can be very useful when you are
	trying to debug the code and want to see if results are
	reproducible.  Sometimes I set up my code so that it asks for a
	seed, but if the seed I give is negative, it calls pid() to get a
	random seed.  That way, I can do many runs without having to think
	up new seeds, and when I want to reproduce runs with a known seed I
	can also do that.

	Call pi1 with the argument 50 and 500, and send the output to a
	file.  Each line will have 50 (the number of darts) and the
	estimate of pi/4 (the value of the integral).

	Hist

	It is interesting to look at the distribution of the estimates.  One
	way to do this is by making a histogram.  You break the x-interval
	up into different ranges and then count the number of times a
	value falls in each range.  hist is a program for doing this.  I
	find it very useful.  You are welcome to copy the source code from
	~sg/src/hist.c.  If you have a file with a bunch of numbers to
	histogram, just type hist.c <file.  By default, there will be 10
	bins.

	Here is an example:
	pi1 50 500 |awk '{print $2}' |hist

	You can adjust the range or width or number of bins.  You can also
	have the output prepared for axis.  Here is another example:

	pi1 50 500 |awk '{print $2}' |hist -n 20 -g |axis |xplot

	Here is a reminder of all the hist options:
	*       option -n specifies number of bins                      *
	*       option -s specifies size of bin                         *
	*       option -x specifies lower and upper limits of histogram *
	*       option -g specifies output in form suitable for "graph" *
	*                       or "axis"                               *
	*                 if -g 1 we get a full bar for each bin        *
	*                 if -g 2 the bars don't go down to the origin  *
	*       option -m specifies the line type for use with "axis"   *


	Aver

	Aver is another useful program you can copy.  For our purposes you
	only have to know that if you have a file with x and y on each
	line, if you type
	aver file
	then for any lines that have the same x-value aver will average all
	the corresponding y-values.  It will print x y_average y_error.
	This can be very useful when you want averages.  It is also why I
	arranged the output of pi1 so that the number of darts was the
	first word of output on each line.

	Write a shell script to put the results pi1 with a number of
	different trial sizes in a file and then use aver to average the
	results and make a plot of the average value of pi vs. size.
	You can also write an awk script to print the absolute value of the
	error.  Here is the awk commaand:
	{ err=3.1415926-$2; if(err<0) err= (-err); print $1,err}
	Now you many apply this command to your file before sending it to
	aver and then plot the error vs size.  Use logarithmic scale for
	the x and y axes.  Can you verify that the error falls like
	1/sqrt(size) ?

	Higher Dimensional Integrals

	Monte Carlo methods like hit or miss and the sample mean methods
	have errors that fall like n^{-1/2}.  This is not very good.  To
	cut the error in half you need to do 4 times a much work.  For
	Simpson's rule the error goes like n^{-4}, so if you do twice a
	much work, the error falls by a factor of 16.  However, for 
	multidimensional integrals, a method with an error like h^a,
	where h is the spacing between grid points, will have an error like
	n^{-a/d} because the number of points on the grid goes like h^{-d}.
	As an example, if a=4 and d=10, a/d=0.4 so the error is falling
	more slowly than n^{-1/2}.  So, for high dimension integrals,
	Monte Carlo methods may be bad, but other methods are worse.