# S functions for inference from iterative simulation # Written 1 July 1991 by Andrew Gelman, Dept. of Statistics, UC Berkeley # Please mail all comments/questions to gelman@stat.berkeley.edu # Software can be freely used for non-commercial purposes and freely # distributed. Submitted to statlib on 29 April 1992. # Derivations appear in the article, "Inference from iterative simulation," # by A. Gelman and D. B. Rubin, in "Statistical Science," to appear. # How to Use These Programs # # Preparation: The results of m multiple sequences, each of length 2n, # of iterative simulation (e.g., the Metropolis algorithm or the # Gibbs sampler) used to simulate a probability distribution. # The starting points of the simulations should be drawn from a # distribution that is overdispersed relative to the "target" # distribution you are trying to simulate. # # At each iteration, all scalar parameters or summaries of interest # should be stored. The result will be an 2n by m matrix for each # scalar summary. # # To run: Use gpar (or gpar.log or gpar.logit) to do output analysis # on the simulations for each scalar summary of interest. For # each scalar summary r, gpar produces two results: # # a. Estimates of the distribution of r, # b. The estimated potential scale reduction; # i.e., the factor by which the estimated t intervals for r # might decrease if the iterative simulation is continued # forever. As n increases, the scale reduction approaches 1. # The following S functions are included: # gpar (Gibbs-parallel) takes a matrix r of simulations: each of # m columns is an independent iteratively simulated sequence of length # 2n. The first n steps of each sequence are discarded. # Output is a list of three vectors: # # post: (2.5%, 50%, 97.5%) quantiles for the target distribution # based on the Student-t distribution # quantiles: (2.5%, 25%, 50%, 75%, 97.5%) empirical quantiles of # the mn simulated values. # confshrink: 50% and 97.5% quantiles of a rough upper bound on # how much the confidence interval of "post" will shrink # if the iterative simulation is continued forever. # # If both components of confshrink are not near 1, you should # probably run the iterative simulation further. # gpar.log and gpar.logit are versions of gpar to be used for scalar # summaries that are all-positive or are constrained to lie # between 0 and 1. ############################################################################# # SUBROUTINES cov _ function (a,b) { m _ length(a) ((mean ((a-mean(a))*(b-mean(b)))) * m)/(m-1)} logit _ function (x) {log(x/(1-x))} invlogit _ function (x) {1/(1+exp(-x))} col.means _ function(mat) { ones _ matrix(1, nrow = 1, ncol = nrow(mat)) ones %*% mat/nrow(mat)} col.vars _ function(mat) { means _ col.means(mat) col.means(mat * mat) - means * means} # Chi-squared degrees of freedom estimated by method of moments # # (Assume A has a gamma sampling distribution and varA is an unbiased # estimate of the sampling variance of A.) chisqdf _ function (A, varA) {2*(A^2/varA)} ############################################################################# # MAIN PROGRAM gpar _ function (r) { alpha _ .05 # 95% intervals m _ ncol(r) x _ r [(nrow(r)/2+1):nrow(r),] # second half of simulated sequences n _ nrow(x) # We compute the following statistics: # # xdot: vector of sequence means # s2: vector of sequence sample variances (dividing by n-1) # W = mean(s2): within MS # B = n*var(xdot): between MS. # muhat = mean(xdot): grand mean; unbiased under strong stationarity # varW = var(s2)/m: estimated sampling var of W # varB = B^2 * 2/(m+1): estimated sampling var of B # covWB = (n/m)*(cov(s2,xdot^2) - 2*muhat*cov(s^2,xdot)): # estimated sampling cov(W,B) # sig2hat = ((n-1)/n))*W + (1/n)*B: estimate of sig2; unbiased under # strong stationarity # quantiles: emipirical quantiles from last half of simulated sequences xdot _ as.vector(col.means(x)) s2 _ as.vector(col.vars(x)) W _ mean(s2) B _ n*var(xdot) muhat _ mean(xdot) varW _ var(s2)/m varB _ B^2 * 2/(m-1) covWB _ (n/m)*(cov(s2,xdot^2) - 2*muhat*cov(s2,xdot)) sig2hat _ ((n-1)*W + B)/n quantiles _ quantile (as.vector(x), probs=c(.025,.25,.5,.75,.975)) if (W > 1.e-8) { # non-degenerate case # Posterior interval post.range combines all uncertainties # in a t interval with center muhat, scale sqrt(postvar), # and postvar.df degrees of freedom. # # postvar = sig2hat + B/(mn): variance for the posterior interval # The B/(mn) term is there because of the # sampling variance of muhat. # varpostvar: estimated sampling variance of postvar postvar _ sig2hat + B/(m*n) varpostvar _ (((n-1)^2)*varW + (1+1/m)^2*varB + 2*(n-1)*(1+1/m)*covWB)/n^2 post.df _ chisqdf (postvar, varpostvar) post.range _ muhat + sqrt(postvar) * qt(1-alpha/2, post.df)*c(-1,0,1) # Estimated potential scale reduction (that would be achieved by # continuing simulations forever) has two components: an estimate and # an approx. 97.5% upper bound. # # confshrink = sqrt(postvar/W), # multiplied by sqrt(df/(df-2)) as an adjustment for the # width of the t-interval with df degrees of freedom. # # postvar/W = (n-1)/n + (1+1/m)(1/n)(B/W); we approximate the sampling dist. # of (B/W) by an F distribution, with degrees of freedom estimated # from the approximate chi-squared sampling dists for B and W. (The # F approximation assumes that the sampling dists of B and W are independent; # if they are positively correlated, the approximation is conservative.) varlo.df _ chisqdf (W, varW) confshrink.range _ sqrt (c(postvar/W, (n-1)/n + (1+1/m)*(1/n)*(B/W) * qf(.975, m-1, varlo.df)) * post.df/(post.df-2)) list(post=post.range, quantiles=quantiles, confshrink=confshrink.range) } else { # degenerate case: all entries in "data matrix" are identical list (post=muhat*c(1,1,1), quantiles=quantiles, confshrink=c(1,1)) } } ############################################################################## gpar.log _ function (r) { gp _ gpar(log(r)) list (post=exp(gp\$post), quantiles=exp(gp\$quantiles), confshrink=gp\$confshrink)} gpar.logit _ function (r) { gp _ gpar(logit(r)) list (post=invlogit(gp\$post), quantiles=invlogit(gp\$quantiles), confshrink=gp\$confshrink)} gpar.vec _ function (a) { # monitors the convergence of the gibbs sampler # a[i,j,k], i=1,...,n, j=1,...,m, k=1,...,d d _ dim(a)[3] check _ NULL for (k in (1:d)){ gp _ gpar(a[,,k]) check _ rbind (check, gp\$confshrink) } check } "do.summary"<- function(a, v = c(0, 0.01, 0.025, 0.05, 0.1, 0.25, 0.5, 0.75, 0.9, 0.95, 0.975, 0.99, 1), coef.names = paste(1:(dim(a)[3]))) { x <- apply(a, 3, quantile, probs = v) dimnames(x) <- list(paste(round(100 * v, 0), "%", sep = ""), coef.names ) t(x) }