November 16, 2004 --- Class 23---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.  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 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 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) ?

	Next homework

	It is convenient use the split utility, or my blockavg program for
	the sample mean problem in the next homework.  We discussed how
	these may be used.

	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.  One useful example where this works in a exponential
	distribution.  We were able to show that 

	P(x) = 1 - exp(-x/lambda)
	
	Setting P(x)=r and solving for r, we find

	x = - lambda * ln(1-r) 

	Since 1-r and r are both uniformly distributed between 0 and 1,
	we can simplify the above and replace ln(1-r) by ln(r).
	Thus, to generate numbers with an exponential distribution, we
	first generate a uniform random variable between 0 and 1.  Then
	
	x = -lambda * ln(r)

	has an exponential distribution.

	In the next class, we will use Mathematica to test this.