Software and examples for combining climate model output

This page provides a tutorial to the Software and data sets used in: Tebaldi, Smith Nychka and Mearns (2004) Quantifying Uncertainty in Projections of Regional Climate Change: A Bayesian Approach to the Analysis of Multimodel Ensembles

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.cgd.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()
 [1] "A2.DJF.Y"    "A2.JJA.Y"    "B2.DJF.Y"    "B2.JJA.Y"    "DJF.X"
 [6] "DJF.X0"      "DJF.lambda0" "JJA.X"       "JJA.X0"      "JJA.lambda0"
[11] "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:
(scenario).season.variable The X and Y variables are matrices where rows index regions and columns index models. Thus A2.DJF.Y is a 22X9 matrix of the future projections for the A2 scenario in winter. Part of the reason that the data are grouped together is that the Gibbs sampler is more efficient sampling for all regions simultaneously. (For those familiar with R programming we have eliminated an explicit loop over regions using matrix multiplications.) X0 and lambda0 refer to the climate and associated precision based on observed data. The regions are coded with three letters e.g. CAS= Central Asia, WNA=Western North America.

An example:

Source the two files in R as indicated above and now in a R session

# 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"])