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