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[]