-
Notifications
You must be signed in to change notification settings - Fork 0
/
determining-convergence.Rmd
138 lines (109 loc) · 6.58 KB
/
determining-convergence.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
---
title: "Determining Convergence"
output: rmarkdown::html_vignette
vignette: >
%\VignetteIndexEntry{Determining Convergence}
%\VignetteEngine{knitr::rmarkdown}
%\VignetteEncoding{UTF-8}
bibliography: references.bib
---
```{r, include = FALSE}
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>",
dev='png',
fig.width=7,
fig.height=5.5 # ,
# dev.args=list(antialias = "none")
)
```
```{r setup}
library(carbondate)
```
# Introduction
A few tools to determine how well the MCMC for the various approaches have converged.
## Examining the Gelman-Rubin convergence diagnostic
This is often used to evaluate MCMC convergence. It compares the _between-chains_ variance with the _within-chains_ variance of the model parameters for multiple MCMC chains. If the MCMC has converged to the target posterior, then these values should be similar.
To assess convergence of our methods, we apply it to the individual posterior calendar age estimates of the ^14^C samples. We have $N$ such samples, each of which has a individual sequence of sampled calendar age values that are stored within the overall MCMC. We calculate the potential scale reduction factor (PSRF) for each individual sample's calendar age sequence, generating $N$ PSRF values. If the MCMC has converged to the target posterior distribution, then each of these PSRFs should be close to 1.
In the first case, the relevant MCMC function can be run multiple times. This generate different chains.
```{r calculate_gr_multiple, results=FALSE, eval=FALSE}
all_outputs <- list()
for (i in 1:3) {
set.seed(i + 1)
all_outputs[[i]] <- PolyaUrnBivarDirichlet(
kerr$c14_age, kerr$c14_sig, intcal20, n_iter = 1e4)
}
PlotGelmanRubinDiagnosticMultiChain(all_outputs)
```
```{r, out.width= "100%", echo = FALSE}
knitr::include_graphics("figures-convergence/calculate_gr_multiple-1.png")
```
It can also be calculated by taking a single MCMC run, and splitting it into multiple parts to compare the _within-segment_ variance with the _between-segment_ variance for each calendar age observation.
```{r calculate_gr, results=FALSE, eval=FALSE}
set.seed(3)
output <- PolyaUrnBivarDirichlet(
kerr$c14_age, kerr$c14_sig, intcal20, n_iter = 2e4)
PlotGelmanRubinDiagnosticSingleChain(output, n_burn = 5e3)
```
```{r, out.width= "100%", echo = FALSE}
knitr::include_graphics("figures-convergence/calculate_gr-1.png")
```
As you can see, even with a few iterations (where we would expect the result not to have converged yet) the PSRF values are close to one.
## Examining the predictive distribution or posterior occurrence rate
When calibrating multiple ^14^C determinations, primary interest will frequently be in the summarised (predictive) calendar age estimate if using the Bayesian non-parametric method, or the posterior occurrence rate if using the Poisson process approach, rather than the age of any individual sample. Information on the predictive estimate is encapsulated in the model parameters relating to the underlying clusters (e.g., weights and distributions). While the occurrence rate is defined by the locations of the changepoints, and the segment heights, in the Poisson process model.
Unfortunately, the chains storing these parameters are not suitable for the Gelman-Rubin diagnostic. In the case of the predictive Bayesian non-parametric estimate, the number and identity of the clusters stored in the MCMC can change with each iteration (as clusters drop-in and out, or are relabelled). In the case of the Poisson process occurrence rate, the number and labelling of changepoints varies throughout the MCMC. We therefore provide a further diagnostic, based upon assessing the predictive calendar age estimate or posterior occurrence rate, which may be a much more useful indicator of MCMC convergence.
### Visually comparing multiple runs
Running the functions a few time with different random number seeds can give an idea of how many iterations are needed for convergence. If the MCMC has converged, then each run should lead to a similar result for the predictive density (or the posterior occurrence rate in the case of the POisson process model). For example,
```{r calculate_polya_kerr, results=FALSE, eval=FALSE}
outputs <- list()
for (i in 1:3) {
set.seed(i+1)
outputs[[i]] <- PolyaUrnBivarDirichlet(
rc_determinations = kerr$c14_age,
rc_sigmas = kerr$c14_sig,
calibration_curve=intcal20,
n_iter = 1e4)
outputs[[i]]$label <- paste("Seed =", i)
}
PlotPredictiveCalendarAgeDensity(
outputs, n_posterior_samples = 500, denscale = 2, interval_width = "1sigma")
```
```{r, out.width= "100%", echo = FALSE}
knitr::include_graphics("figures-convergence/calculate_polya_kerr-1.png")
```
As you can see, in this case (255 determinations collated by Kerr and McCormick [@kerr2014]) the different runs do not have similar outputs, so more iterations would be needed to ensure convergence.
In contrast, if we run a much simpler example (that of artificial data comprised of two normals), we can see that convergence appears to be achieved in a small number of iterations.
```{r calculate_polya_normals, results=FALSE, eval=FALSE}
outputs <- list()
for (i in 1:3) {
set.seed(i + 1)
outputs[[i]] <- PolyaUrnBivarDirichlet(
rc_determinations = two_normals$c14_age,
rc_sigmas = two_normals$c14_sig,
calibration_curve=intcal20,
n_iter = 1e4)
outputs[[i]]$label <- paste("Seed =", i)
}
PlotPredictiveCalendarAgeDensity(
outputs, n_posterior_samples = 500, denscale = 2, interval_width = "1sigma")
```
```{r, out.width= "100%", echo = FALSE}
knitr::include_graphics("figures-convergence/calculate_polya_normals-1.png")
```
This approach to assessing convergence can be taken with either the Bayesian non-parametric method, or the Poisson process modelling.
### Examining the Kullback--Leibler divergence (Bayesian non-parametrics only)
We also provide a further diagnostic specifically for the Bayesian non-parametric approach (it is not applicable for the Poisson process model). This diagnostic gives a measure of the difference between an initial (baseline) predictive density and the predictive density as the MCMC progresses.
```{r calculate_kld, results=FALSE, eval=FALSE}
set.seed(50)
output <- WalkerBivarDirichlet(
rc_determinations = kerr$c14_age,
rc_sigmas = kerr$c14_sig,
calibration_curve=intcal20,
n_iter = 1e5)
PlotConvergenceData(output)
```
```{r, out.width= "100%", echo = FALSE}
knitr::include_graphics("figures-convergence/calculate_kld-1.png")
```
It can give an idea of convergence as well as which iteration number to use for `n_burn` when calculating the predictive density (by default set to half the chain).
# References