Skip to content

Commit

Permalink
Merge pull request #46 from TJHeaton/development
Browse files Browse the repository at this point in the history
Acceptance by CRAN as v1.0.1
  • Loading branch information
TJHeaton committed Jan 31, 2024
2 parents 7bb32f4 + a057006 commit 070f528
Show file tree
Hide file tree
Showing 15 changed files with 36 additions and 15 deletions.
7 changes: 3 additions & 4 deletions CRAN-SUBMISSION
Original file line number Diff line number Diff line change
@@ -1,4 +1,3 @@
Version: 1.0.0
Date: 2024-01-25 12:00:37 UTC
SHA:
61ef58cf978a4b4ef7dc3f2cb66a4dc5ee7039e8
Version: 1.0.1
Date: 2024-01-29 15:41:13 UTC
SHA: b8e397982597cdc1eef472752ecde7af1e5f250d
2 changes: 1 addition & 1 deletion DESCRIPTION
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
Package: carbondate
Title: Calibration and Summarisation of Radiocarbon Dates
Version: 1.0.0
Version: 1.0.1
Authors@R: c(
person("Timothy J", "Heaton", , "T.Heaton@leeds.ac.uk", role = c("aut", "cre", "cph"),
comment = c(ORCID = "0000-0002-9994-142X")),
Expand Down
4 changes: 4 additions & 0 deletions NEWS.md
Original file line number Diff line number Diff line change
@@ -1,3 +1,7 @@
# carbondate 1.0.1

* Fixing CRAN comments

# carbondate 1.0.0

* First release
Expand Down
2 changes: 1 addition & 1 deletion README.Rmd
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,7 @@ knitr::opts_chunk$set(

<!-- badges: start -->
[![R-CMD-check](https://github.com/TJHeaton/carbondate/actions/workflows/R-CMD-check.yaml/badge.svg)](https://github.com/TJHeaton/carbondate/actions/workflows/R-CMD-check.yaml)
![CRAN status](https://www.r-pkg.org/badges/version/carbondate)](https://CRAN.R-project.org/package=carbondate)
[![CRAN status](https://www.r-pkg.org/badges/version/carbondate)](https://CRAN.R-project.org/package=carbondate)
<!-- badges: end -->

An R package to analyse, and summarise, multiple radiocarbon (^14^C) determinations. The package provides two linked (but distinct) Bayesian approaches that can both be used to obtain rigorous and robust alternatives to summed probability distributions (SPDs):
Expand Down
4 changes: 2 additions & 2 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -6,8 +6,8 @@
<!-- badges: start -->

[![R-CMD-check](https://github.com/TJHeaton/carbondate/actions/workflows/R-CMD-check.yaml/badge.svg)](https://github.com/TJHeaton/carbondate/actions/workflows/R-CMD-check.yaml)
![CRAN
status](https://www.r-pkg.org/badges/version/carbondate)\](<https://CRAN.R-project.org/package=carbondate>)
[![CRAN
status](https://www.r-pkg.org/badges/version/carbondate)](https://CRAN.R-project.org/package=carbondate)
<!-- badges: end -->

An R package to analyse, and summarise, multiple radiocarbon
Expand Down
4 changes: 3 additions & 1 deletion cran-comments.md
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,9 @@ We have also:

* Altered documentation of functions (and vignettes) to incorporate above plotting changes (and fix a few typos to ensure parameter naming consistency across outputs).

* Fixed typos in vignettes and run MCMC in "Non-parametric-summed-density" for longer to show converged result
* Fixed typos in vignettes

* Run MCMC in "Non-parametric-summed-density" for longer to show converged result, shortened other vignettes to get under 10min runtime

* Switched antialiaising back to default on png on vignettes.

Expand Down
Binary file modified man/figures/README-plot_density-1.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file modified man/figures/README-plot_mean_rate-1.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
3 changes: 2 additions & 1 deletion vignettes/Against_SPDs.Rmd
Original file line number Diff line number Diff line change
Expand Up @@ -118,6 +118,7 @@ polya_urn_output <- PolyaUrnBivarDirichlet(
rc_determinations = two_normals$c14_age,
rc_sigmas = two_normals$c14_sig,
calibration_curve=intcal20,
n_iter = 1e4,
show_progress = FALSE)
two_normals_DPMM <- PlotPredictiveCalendarAgeDensity(
Expand All @@ -130,7 +131,7 @@ par(new = TRUE,
yaxs = "i",
mar = c(5, 4.5, 4, 2) + 0.1,
las = 1)
xlim <- rev(range(two_normals_DPMM[[1]]$calendar_age))
xlim <- rev(range(two_normals_DPMM[[1]]$calendar_age_BP))
ylim <- c(0, 3 * max(two_normals_DPMM[[1]]$density_mean))
plot(
NULL,
Expand Down
25 changes: 20 additions & 5 deletions vignettes/determining-convergence.Rmd
Original file line number Diff line number Diff line change
Expand Up @@ -34,7 +34,7 @@ To assess convergence of our methods, we apply it to the individual posterior ca

In the first case, the relevant MCMC function can be run multiple times. This generate different chains.

```{r calculate_gr_multiple, results=FALSE}
```{r calculate_gr_multiple, results=FALSE, eval=FALSE}
all_outputs <- list()
for (i in 1:3) {
set.seed(i + 1)
Expand All @@ -43,16 +43,22 @@ for (i in 1:3) {
}
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}
```{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.

Expand All @@ -66,7 +72,7 @@ Unfortunately, the chains storing these parameters are not suitable for the Gelm

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}
```{r calculate_polya_kerr, results=FALSE, eval=FALSE}
outputs <- list()
for (i in 1:3) {
set.seed(i+1)
Expand All @@ -80,12 +86,15 @@ for (i in 1:3) {
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}
```{r calculate_polya_normals, results=FALSE, eval=FALSE}
outputs <- list()
for (i in 1:3) {
set.seed(i + 1)
Expand All @@ -99,14 +108,17 @@ for (i in 1:3) {
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}
```{r calculate_kld, results=FALSE, eval=FALSE}
set.seed(50)
output <- WalkerBivarDirichlet(
rc_determinations = kerr$c14_age,
Expand All @@ -116,6 +128,9 @@ output <- WalkerBivarDirichlet(
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).

Expand Down
Binary file added vignettes/figures-convergence/calculate_gr-1.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.

0 comments on commit 070f528

Please sign in to comment.