Produces plots for output and diagnostics from MCMC simulations (D.A. Murray).
Options
PRINT = string tokens |
Controls printed output (summary); default * |
|---|---|
PLOT = string tokens |
Controls the type of plot (trace, density, autocorrelation, gelmanrubin); default trac |
ARRANGEMENT = string tokens |
Specifies whether to draw the plots individually or 4 to a page (single, multiple); default sing |
START = scalar |
Start iteration number for plots |
END = scalar |
End iteration number for plots |
MAXLAG = scalar |
Maximum lag for autocorrelation plots; default 50 |
BANDWIDTH = scalar |
The bandwidth value to be used for the density plots. |
GRMETHOD = scalar |
Controls the method of the Gelman-Rubin diagnostic plot (gr, bgr); default bgr |
BINWIDTH = scalar |
Number of values in each bin in the Gelman-Rubin plot; default 50 |
USEALLSAMPLES = text |
Whether to use all the samples for Gelman-Rubin plot, or to discard the first half of the observations (yes, no); default no |
Parameter
SIMULATIONS = pointers |
List of pointers containing simulations, one for each Markov chain |
|---|
Description
BGPLOT produces plots for output and diagnostics from Markov Chain Monte Carlo (MCMC) simulations. The procedure can be used after running a Bayesian MCMC analysis using BGXGENSTAT. Alternatively, data from CODA files produced by WinBUGS that are imported using BGIMPORT can be displayed using BGPLOT.
The data are supplied using the SIMULATIONS parameter as list of pointers, where each pointer contains the simulations and associated information for a Markov chain. The data for BGPLOT can be taken directly from the BGXGENSTAT or BGIMPORT procedures. The SIMULATIONS parameter corresponds to the SIMULATIONS parameters of BGXGENSTAT and BGIMPORT.
BGPLOT produces four types of plot which can be selected using the PLOT option. The ARRANGEMENT option controls whether the plots are each drawn on separate pages or four to a page in a two-by-two arrangement.
The trace setting of PLOT produces a trace plot for each monitored node, where every chain for the monitored nodes is superimposed on the same plot.
The density setting produces a kernel density smooth for each monitored node. The bandwidth for the density plot can be supplied using the BANDWIDTH option.
The autocorrelation setting generates an autocorrelation plot for each monitored node, where every chain for the monitored nodes is superimposed on the same plot. You can specify the maximum lag for which the autocorrelation is calculated using the MAXLAG option (default 50).
The gelmanrubin setting produces the Gelman-Rubin “Potential Scale Reduction Factor” (PSRF) diagnostic. The PSRF compares the between and within variances of multiple chains. BGPLOT produces a plot of the PSRF, and a second plot of the between and within variances. Convergence is assumed to occur when the PSRF is close to one, and the between and within variances stabilise around the same value. You can display the convergence diagnostic by setting option PRINT=summary. The method used for the diagnostic can be controlled using the GRMETHOD option. The gr setting uses the original Gelman-Rubin PSRF diagnostic, while the bgr setting uses the Gelman-Rubin-Brooks version, which interprets the diagnostic as a ratio of interval lengths rather than a variance ratio. The BINWIDTH option specifies the number of values in each bin of the Gelman-Rubin plot (default 50). Usually the first half of the observations is discarded in the Gelman-Rubin plot, but you can set option USEALLSAMPLES=yes to use them all.
Options: PRINT, PLOT, ARRANGEMENT, START, END, MAXLAG, BANDWIDTH, GRMETHOD, BINWIDTH, USEALLSAMPLES.
Parameter: SIMULATIONS.
Method
The autocorrelations are calculated using the CORRELATE directive and the densities are evaluated using the KERNELDENSITY procedure. Gelman & Rubin’s diagnostic (1992) monitors convergence of MCMC output from parallel chains which each have different starting points that are overdispersed with respect to the target distribution. The details for the calculation of the potential scalar reduction factor (R) is described in Gelman & Rubin (1992). The Gelman-Rubin-Brooks diagnostic is an alternative method to the R statistic devised by Gelman & Rubin where R is calculated by the ratio of the central 80% width of the pooled intervals by the average central 80% width of the within intervals, see Brooks & Gelman (1998).
Action with RESTRICT
Any data restrictions will be ignored.
References
Brooks, S.P. & Gelman, A. (1998). General Methods for Monitoring Convergence of Iterative Simulations. Journal of Computational and Graphical Statistics, 7, 434-455.
Gelman, A. & Rubin, D. (1992). Inference from Iterative Simulation Using Multiple Sequences. Statistical Science, 7, 457-511.
See also
Procedures: BGIMPORT, BGXGENSTAT.
Commands for: Bayesian methods.
Example
CAPTION 'BGPLOT Example','Example from WinBUGS.',\
!t('Binary dose-response data published by Bliss (1935), in which',\
'the numbers of beetles killed after 5 hour exposure to carbon',\
'disulphide at N=8 different concentrations are recorded.',\
'Data analysed by logit model using WinBUGs.'),\
!t('Note that WinBUGS or OpenBUGs must be installed on your system',\
'to run this example.'); style=meta,3(plain)
" Specify WinBUGS model "
TEXT tmodel; VALUES=!t(\
'model',\
'{',\
' for( i in 1 : N ) {',\
' r[i] ~ dbin(p[i],n[i])',\
' logit(p[i]) <- alpha + beta * (x[i] - mean(x[]))',\
' rhat[i] <- n[i] * p[i]',\
' phat[i] <- r[i]/n[i]',\
' }',\
' beta ~ dnorm(0.0,0.00001)',\
' alpha ~ dnorm(0.0,0.00001)',\
'}')
" Write the WinBUGS model to a file "
OPEN 'beetles-model.txt'; CHAN=2; FILE=output
PRINT [CHANNEL=2; IPRINT=*; SQUASH=yes] tmodel; JUST=left
CLOSE 2; FILETYPE=output
" Data "
SCALAR N; VALUE=8
VARIATE x; VALUES=!(1.6907,1.7242,1.7552,1.7842,1.8113,1.8369,1.8610,1.8839)
VARIATE n; VALUES=!(59,60,62,56,63,59,62,60)
VARIATE r; VALUES=!(6,13,18,28,52,53,61,60)
TEXT dnames; VALUES=!t(N,x,n,r)
POINTER [VALUES=N,x,n,r] data
" Initial values "
TEXT ipnames; VALUES=!T('alpha','beta')
POINTER [NVALUES=ipnames] init
SCALAR init1[1...2]; VALUE=100,-100
SCALAR init2[1...2]; VALUE=-100,100
" Run WinBugs"
BGXGENSTAT [PRINT=node; MONITOR=ipnames;\
MODELFILE='beetles-model.txt';\
DATA=data; IDATANAMES=dnames; NSAMPLES=5000;\
INAMES=ipnames; NCHAIN=2] init1,init2; SIMULATIONS=sim[1...2]
" Plots "
BGPLOT [PRINT=summary; PLOT=trace,density,auto,gelm;\
ARRANGEMENT=multiple] sim[]