November 7, 2006 --- Class 20---Hit or Miss Method, Inverse Transform Method

Activities:

	Hit-or-Miss Method

	The Hit-or-Miss method and Sample Mean methods are described in 
	Section 11.2.  We reviewed the hit-or-miss method and detailed the
	sample mean method.  We'll try the first method and see how well it does.
	However, first I showed an example of the converge of the trapezoid
	rule and Simpson's run.

	Trapezoid rule and Simpson's rule

	I did this example on my work computer and moved all the work to
	Nations using a tar (tape archive) file.  To create a tar file, use

	tar -cvf myarchive.tar  files

	or

	tar -cvf myarchive.tar  directory

	c is for create, v is for verbose and f indicated that the next
	argument is the name of the tape archive file to be created.  Then
	comes a list of files and/or directories that are to be include in the
	archive.  To extract the files from the archive, you may use the
	command:

	tar -xvf myarchive.tar

	This unpacks the files from myarchive.tar and puts them in the current
	directory.  If you used the second form above, then a directory will
	be created and the files placed there.  
	
	In ~sg/integration_test/double_precision, you will find source code 
	and output files.  The source code has been modified so we can see
	how the method converges.  The integral of cos(x) from x=0 to x=pi is
	calculated by xqsimp2.  Here is the output:

	1       7.853981633974e-01      3.333333333333e+29
	2       9.480594489685e-01      1.002279877492e+00
	3       9.871158009728e-01      1.000134584974e+00
	4       9.967851718862e-01      1.000008295524e+00
	5       9.991966804851e-01      1.000000516685e+00
	6       9.997991943200e-01      1.000000032265e+00
	7       9.999498000921e-01      1.000000002016e+00
	8       9.999874501175e-01      1.000000000126e+00
	9       9.999968625353e-01      1.000000000008e+00
	10      9.999992156342e-01      1.000000000000e+00

	The first column is the level at which we call the routine, i.e., with
	2^j points in the integration range.  The second column is the
	trapzoid rule and the third column is for Simpson's rule.  Note that
	the j=1 value for Simpson's rule is bogus since we need to go at least
	to j=2 to use Simpson's rule.  In cos.ax, you will find the error
	plotted vs. j.  The crosses are for the trapezoid rule, the diamonds
	are for Simpson's rule.

	A second example integrates sqrt(1-x^2) from 0 to 1.  Curiously, it
	looks very different.  This is because of the singularity in the
	derivative of the integrand at x=1.  Simpson's rule does not converge
	with a higher power of j.  Plot qsimp.ax, or plot both graphs
	together.

	Testing Hit-or-Miss Method

	Program ~sg/src/pi3.c implements the hit-or-miss method using a
	feedback shift register random number generator.
	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 pi3 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:
	pi3 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:

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

	Here are 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 pi3 so that the number of darts was the
	first word of output on each line.

	Write a shell script to put the results pi3 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) ?

	Alternatively, you can get the variance by multiplying the last column
	by the square root of the number of batches you averaged over in each
	case.  I showed a log log graph that verified tha the error falls like
	the square root of the number of trials.  

	Inverse Transform Method

	Sometimes we are interested in generating random numbers with a non
	uniform probabilty density.  There are a number of ways of doing
	this and we shall study a few.  This is covered in CSM section
	11.5, starting on page 358.  The first method we will study is
	called inverse transform.  You may be used to the idea of the
	probability density function.  For a real variable, we have
	$\int_{-\Infinity}^{+\Infinity} p(x) dx = 1 $.

	We can define the monotonically increasing function called the
	cumulative probability distribution by
	$ P(x) = \int_{-\Infinity}^{x} p(x') dx' $.

	Of course, $ d P(x)/dx = p(x)$.  P(x) is a number between 0 and 1.
	If we can solve the equation $ P(x)=r$  for $x(r)$, then if we
	generate uniform random numbers r, then the numbers x(r) have the
	distribution we want.  So, the key is we start with p(x), and
	the method only is useful if we can find P(x) and then solve
	P(x)=r for x.

	In our next class, we will try some examples.