This function performs the Markov Chain Monte Carlo sampling for the models from the paper Tebaldi et al. (2004). The default values for the Gibbs sampler are the same as those used in the analysis of the paper. For a climate scientist willing to learn R, this will allow you to investigate MCMC applications in the privacy of your office!
REA.Gibbs(X, X0, lambda0, Y, N = 5000, save.every = 50, burn.in = 100 * save.every)
X 
A matrix of climate model results for present climate. Rows index regions and columns index models. 
X0 
A vector or one column matrix of observed climate for the regions. 
lambda0 
A vector or one column matrix of the precisions ( reciprocal variances for observed climate for the regions. 
Y 
A matrix of climate model results for future climate. Rows index regions and columns index models. 
N 
The size of the returned sample 
save.every 
The spacing in the chain between members of the sample 
burn.in 
The number of initial members in the chain to skip before starting to sample. 
The function implements a Gibbs sampler for the Gaussianbased statistical model for present and future model biases as described in Quantifying Uncertainty in Projections of Regional Climate Change: A Bayesian Approach to the Analysis of Multimodel Ensembles (2004), Tebaldi, Smith Nychka and Mearns. http://www.cgd.ucar.edu/stats/pub/nychka/manuscripts/tebaldiREA.pdf Technical details are given in www.cgd.ucar.edu/stats/pub/nychka/manuscripts/tebaldisupp.pdf
The idea is to generate a very long multivariate "time series". The components of the time series are the parameters in the statistical model and the stationary distribution of this series is the Bayesian posterior that we want. The discrete time step in this case is just one full iteration of the Gibbs sampler and the total length of the series is the NGIBBS variable within the function.
Note that the exact posterior density function is never computed. The best we can do is to provide a sample that is approximately a random sample from this distribution. The density can then be estimated by a histogram from this sample. Moreover the densities of any bivariate, multivariate or transformation of the parameters can be estimated in a similar manner. The good news is that one can generate as large a sample as one needs and is only limited by computing time. For the paper we consider samples of size 5000, Certainly adequate for examining a univariate density function. The example below in order to run quickly sets the sample size to a modest 250 (N=250). These samples are accumulated by taking widely spaced members of the original time series. The frequency depends on the particular problem and data but we have an interval of 50 (save.every=50). yields parameter vectors that are much less correlated. Thus one can treat the sample returned form this analysis as being independent draws. Another detail of the sampling is that a burn in period is required for the Gibbs time series to move away from transient effects from particular initial conditions. Again the burn in period depends on the problem and data but we have found a period of 5000 Gibbs iterates to be a conservative choice. In terms of arguments to this function this would correspond to ignoring the first 100 samples if they are 50 iterations apart (50*100=5000). An example below illustrates these features.
mu 
Matrix of sampled values for mu, true value for present climate. Columns index regions and the number of rows is N. 
nu 
Matrix of sampled values for nu, true value for future climate. Columns index regions number of rows is N. 
theta 
Matrix of sampled values for theta the inflation/deflation of precision between present future model biases. Columns index regions the the number of rows is N. 
lambda 
Sampled values for lambda, the model bias precisions.
(precision is the reciprocal of variance) This is array with the indexing
lambda[models, regions,Gibbs samples]. The

Doug Nychka , Claudia Tebaldi
Quantifying Uncertainty in Projections of Regional Climate Change: A Bayesian Approach to the Analysis of Multimodel Ensembles Claudia Tebaldi, Richard L. Smith,Doug Nychka, and Linda Mearns (2003) (in and can be downloaded in PDF format from http://www.cgd.ucar.edu/stats/pub/nychka/manuscripts/tebaldiREA.pdf and the technical supplement http://www.cgd.ucar.edu/stats/pub/nychka/manuscripts/tebaldisupp.pdf
# 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 along with the model output DT< A2.DJF.Y["WNA",]DJF.X["WNA",] hist( look$nu[,"WNA"]look$mu[,"WNA"], xlim=range( DT), xlab="Delta T", freq=FALSE) # adding model output along axis text( DT,rep( .08,9), dimnames(A2.DJF.Y)[[2]], adj=0,srt=90, col=4) points( DT,rep( 0,9), pch="o", cex=2) # I am not even going to try to explain all these arguments! # add a smooth curve estimate of the density  this is to be preferred over # histogram density(look$nu[,"WNA"]look$mu[,"WNA"])> hold lines( hold, col=2, lwd=2) # examining climate change for two adjacent regions: DT1< look$nu[,"WNA"]look$mu[,"WNA"] DT2< look$nu[,"CNA"]look$mu[,"CNA"] title("Estimate of Bivariate temperature change") plot( DT1,DT2, xlab="Western NA", ylab="Central NA") # Here is an illustration of the properties of the raw MCMC time series. # set arguments to return a time series that is 1000 iterate of the #Gibbs sampler  not every 50. Also we are not throwing any of the # beginning. REA.Gibbs( DJF.X, DJF.X0, DJF.lambda0, A2.DJF.Y, N=2000,save.every=1, burn.in=1)>look # the Gibbs time series: matplot(1:2000, cbind(look$mu[,"WNA"], look$nu[,"WNA"]) , type="l") # The few samples we would get with a spacing of 50 seq(1,2000,50)>tt matpoints(tt, cbind(look$mu[tt,"WNA"], look$nu[tt,"WNA"]) , pch="O",col=4, cex=2)