November 16, 2006---Class 22---Metropolis Method, Autocorrelations

Activities:

	Metropolis Method
	-----------------

	Another interesting method was developed by Metropolis, Rosenbluth,
	Rosenbluth, Teller and Teller.  It can be used to generate a random
	walk that has the desired probability distibution.  The key idea is
	that we take the current value of x, take a random trial step, and then
	test the new point's probability.  If it is more probable, we accept
	the new point.  If it less probable, we accept it with probability
	p(x_trial)/p(x).  In the next class a C code will be shown to the 
	class that accomplishes this.  

	The method is discussed in CSM section 11.7, starting on page 435.
	
	The Metropolis method is regarded as one of the 10 greatest numerical
	methods developed in the 20th century.  It's use is quite widespread.
	It is not a particularly good method for generating gaussian random
	numbers, but we can learn a lot about the method, its efficiency and
	data analysis by studing the method in this very simple context.

        I displayed the code to implement the metropolis method and showed how
        detailed balance is achieved by the algorithm.  Because this is such
	an important method, it is important for me that you implement it
	yourself on the next homework assignment.

        Let's explore the efficiency of this method which often
        goes by the name Monte Carlo because of the prominence of random
        numbers.

	Metropolis Optimization
	-----------------------

	Time history, histogram

	Create some numbers with the command
	metropolis 1 5000 0.2 >data
	This should create a run with 5000 numbers of a Gaussian with width
	1.  A maximum step size of +- 0.2 is used.

	Create the time history with the command
	axis a <data |xplot

	Make a histogram with the command
	hist -n 25 -g <data |axis |xplot

	You will notice the time history shows significant correlation from
	one number to the next.  This should not be surprising since the
	maximum difference between x_n and x_{n+1} = 0.2.

	Autocorr

	We can quantify the notion of time correlation by defining a
	quantity called the autocorrelation.  This is defined in CSM on
	page 367 in Eq. (11.52).  The integer j in that equation is called
	the lag.  It is the time difference of the two variables being
	discussed.  In many types of simulations for long lags the
	autocorrelation falls of exponentially.  The characteristic time
	of this falloff is called the autocorrelation time.

	I have a program called autocorr that can compute the
	autocorrelation.  It takes a list of numbers from the standard
	input.  It takes one argument, the maximum lag for which it should
	calculate the autocorrelation.  You can find the source for
	autocorr in ~sg/src.  You will also need statistics.c.  There is a
	Makefile that shows how to make autocorr.

	Example:
	autocorr 20 <data

	will print to the screen the autocorrelation to a maximum lag of
	20.  You may also try
	autocorr 20 <data |axis y l |xplot

	You need to learn something about the error in the
	autocorrelation.  Make a long run of 50000 points.  Put the data in
	a file that I shall call bigrun.

	split
	-----

	A useful utility called split can help with the next exercise.
	If you type
	split -5000  bigrun
	you will find 10 output files, each with 5000 lines.  The will be
	called xaa, xab, xac ....

	You can now run autocorr on each file with the following
	interactive commands:

	foreach i (xa*)
	autocorr 20 <$i >>! autocorr.dat
	end

	The shell will give you continuation prompts for the second and
	third lines since you have started a foreach loop.  The >>! will
	make sure the shell does not complain about appending to a file
	that does not yet exist, if (like me) you like to set the noclobber
	variable.  Your file autocorr.dat will now have 210 lines in it from
	each of the ten runs of the autocorr program.

	If you type
	axis y l x 0 19 <autocorr.dat |xplot
	you will see the results from each run.  I set the maximum x to 19
	so there will not be those pesky straight lines connecting the last
	from from each file with the first point from the next file.

	If you type
	aver autocorr.dat |axis y l e |xplot
	you should see a nice graph with error bars.  This is a fine use of
	the aver program.

	Some more shell scripts

	In ~sg/metropolis you will find some shell and awk scripts that you 
	can use to go from plotting the autocorrelation as a function of
	the lag for fixed delta in the gaussian random walk, to plotting
	for fixed lag, the autocorrelation as a function of delta.  This is
	useful for finding the optimal delta as you are requested to do for
	the homework.  The way these scripts work was discussed in class.