The form of full conditional distributions and the steps of the Gibbs sampler are presented in the Supplement PDF 100K
The Markov Chain Monte Carlo was implemented in the R statistical language. This is a Matlab-like, high level langauage well sutied to statistics and probability. It is freely available from Comprehensive R Archive Network for a variety of platforms. (In particular our calculations were run on a LINUX PC.) We are enthusiatic about R as a mode for data analysis and dissmenating statistical methods (well it is always good to be honest ...). In what follows we will assume the reader is familiar with R. If you are not ,there are excellent tutorials for getting started, see the contributed section under Documents at CRAN. To reproduce the analysis and to peruse the programs transfer the two files REA.data.r REA.Gibbs.r from this directory www.image.ucar.edu/~nychka/REA
We assume that you have the R package. Within an R session:
source("REA.data.temperatures.r") source("REA.Gibbs.r")You need to do this only once provided you save your work space when quitting R. As a a check when you list the contents of the workspace you should see the following.
> ls()  "A2.DJF.Y" "A2.JJA.Y" "B2.DJF.Y" "B2.JJA.Y" "DJF.X"  "DJF.X0" "DJF.lambda0" "JJA.X" "JJA.X0" "JJA.lambda0"  "REA.Gibbs"REA.Gibbs (help file) is the R function for doing the Gibbs sampling and includes comments within the code. The help file gives several examples analyses and one example of is given below. The naming convention for data tables is:
# Run the Gibbs sampler for the A2 future scenario, winter temperatures. # Sample the chain 250 times # The default choices here are to use every 50 th iteration # and "burn" 100 samples in the beginning to allow the chain # to converge to a stationary state. See the help file to change these # options. # # REA.Gibbs( DJF.X, DJF.X0, DJF.lambda0, A2.DJF.Y, N=250)->look # the function as it is running reports the number of samples generated # histogram for posterior for future climate for Western North America hist( look$mu[,"WNA"]) # histogram for posterior for difference in temperatures "delta T" for # Western North America hist( look$nu[,"WNA"]-look$nu[,"WNA"])