## ----include = FALSE---------------------------------------------------------- knitr::opts_chunk$set( collapse = TRUE, comment = "#>", dev='png', fig.width=7, fig.height=5.5 # , # dev.args=list(antialias = "none") ) ## ----setup-------------------------------------------------------------------- library(carbondate) ## ----illustrate_PP, fig.cap = paste("_Modelling sample occurence using a Poisson process with a latent occurence/activity rate. Left Panel: An illustration of a Poisson process with two changepoints in the occurrence rate. The samples (shown as a rug in purple) occur at random calendar times, but proportional to the underlying occurrence rate. There are therefore more samples between (2700, 2300) cal yr BP than other times. Right Panel: Each sample has a ^14^C age (shown as black ticks on the y-axis). We wish to reconstruct the underlying Poisson process rate given only these ^14^C determination._"), out.width="100%", fig.width=10, fig.height=4, echo = FALSE---- oldpar <- par(no.readonly = TRUE) par( mfrow = c(1,2), mgp = c(3, 0.7, 0), xaxs = "i", yaxs = "i", mar = c(5, 4.5, 4, 2) + 0.1, las = 1) set.seed(16) # Artificial start and end dates mincal <- 1950 maxcal <- 3100 t <- seq(mincal, maxcal, by = 1) # Create a simulated lambda (piecewise constant) lambda <- c(rep(0.1, 350), rep(0.45, 400), rep(0.2, length(t) - (400 + 350))) # Choose beta(mutiplier) beta <- 0.2 truerate <- beta * lambda # Sample calendar ages (according to Poisson process) nsamp <- rpois(1, sum(truerate)) truetheta <- sample(t, nsamp, replace = TRUE, prob = truerate) # Create plot showing rate and calendar ages plot(c(min(t), t, max(t)), c(0,lambda,0), lty = 1, lwd = 2, col = "red", type = "l", ylim = c(0, max(lambda) + 0.1), xlim = c(maxcal, mincal) + c(100, -100), yaxs = "i", xlab = "Calendar Age (cal yr BP)", ylab = expression( paste("Occurrence Rate ", lambda, "(", theta, ")" )), main = "Occurence of Samples") rug(truetheta, side = 1, col = "purple", lwd = 2, ticksize = 0.05, quiet = TRUE) abline(v = t[350], lty = 2, lwd = 2) abline(v = t[350 + 400], lty = 2, lwd = 2) # Sample some 14C determinations rc_sigmas <- rep(25, nsamp) rc_determinations <- rnorm(nsamp, mean = approx(intcal20$calendar_age_BP, intcal20$c14_age, truetheta)$y, sd = sqrt(rc_sigmas^2 + (approx(intcal20$calendar_age_BP, intcal20$c14_sig, truetheta)$y)^2) ) # Plot 14C determinations against IntCal20 plot(intcal20$calendar_age_BP, intcal20$c14_age, col = "blue", ylim = range(rc_determinations) + c(-4,4) * max(rc_sigmas), xlim = c(maxcal, mincal) + c(100, -100), xlab = "Calendar Age (cal yr BP)", ylab = expression(paste("Radiocarbon age (", ""^14, "C ", "yr BP)")), type = "l", main = expression(paste(""^14,"C Changepoint Analysis"))) lines(intcal20$calendar_age_BP, intcal20$c14_age + 2 * intcal20$c14_sig, lty = 2, col = "blue" ) lines(intcal20$calendar_age_BP, intcal20$c14_age - 2 * intcal20$c14_sig, lty = 2, col = "blue" ) polygon(c(rev(intcal20$calendar_age_BP), intcal20$calendar_age_BP), c(rev(intcal20$c14_age - 2 * intcal20$c14_sig), intcal20$c14_age + 2 * intcal20$c14_sig), col=rgb(0,0,1,.3), border=NA) rug(rc_determinations, side = 2, quiet = TRUE) # Plot the true rate along the bottom par(new = TRUE) plot(c(min(t), t, max(t)), c(0,truerate,0), lty = 1, col = "red", type = "l", ylim = c(0, 10 * max(truerate)), xlim = c(maxcal, mincal) + c(100, -100), axes = FALSE, xlab = NA, ylab = NA, yaxs = "i") # Reset plotting parameters par(oldpar) ## ----example, out.width="100%"------------------------------------------------ set.seed(15) # Set initial values n_observed <- 40 rc_sigmas <- rep(15, n_observed) # Create artificial rc_determinations calendar_age_range <- c(300, 700) observed_age_range <- c(500, 550) true_theta <- seq( from = observed_age_range[1], to = observed_age_range[2], length = n_observed) intcal_mean <- approx( x = intcal20$calendar_age_BP, y = intcal20$c14_age, xout = true_theta)$y intcal_sd <- approx( x = intcal20$calendar_age_BP, y = intcal20$c14_sig, xout = true_theta)$y rc_determinations <- rnorm( n = n_observed, mean = intcal_mean, sd = sqrt(rc_sigmas^2 + intcal_sd^2)) # Fit the model PP_fit_output <- PPcalibrate( rc_determinations = rc_determinations, rc_sigmas = rc_sigmas, calibration_curve = intcal20, calendar_age_range = calendar_age_range, calendar_grid_resolution = 1, n_iter = 1e5, n_thin = 10, show_progress = FALSE) ## ----meanrate, out.width="100%"----------------------------------------------- # Run plotting function and assign output to PP_posterior_mean_rate_2sigma PP_posterior_mean_rate_2sigma <- PlotPosteriorMeanRate(PP_fit_output) # Note: You can change the PP line widths using the plot_lwd argument # Look at PP_posterior_mean_rate_2sigma (posterior mean and 2 sigma intervals) head(PP_posterior_mean_rate_2sigma) ## ----no_assign_meanrate, eval = FALSE----------------------------------------- # # Run plotting function # PlotPosteriorMeanRate(PP_fit_output) ## ----assign_meanrate, out.width="100%", fig.show='hide'----------------------- # Run plotting function and assign output to PP_posterior_mean_rate_2sigma PP_posterior_mean_rate_2sigma <- PlotPosteriorMeanRate(PP_fit_output) # Will also recreate plot above (but not shown in vignette) # Look at PP_posterior_mean_rate_2sigma (posterior mean and 2 sigma intervals) head(PP_posterior_mean_rate_2sigma) ## ----changepoint_number, out.width="100%"------------------------------------- PlotNumberOfInternalChanges(PP_fit_output) ## ----changepoint_locations, out.width="100%"---------------------------------- PlotPosteriorChangePoints(PP_fit_output) # Can add an n_changes argument, e.g., n_changes = c(2, 3, 4) # if want to condition on a different number of changes ## ----changepoint_rates, out.width="100%"-------------------------------------- PlotPosteriorHeights(PP_fit_output) # As above can add an n_changes argument, e.g., n_changes = c(2, 3, 4) # if want to condition on a different number of changes ## ----PP_plot_individual, out.width="100%"------------------------------------- PlotCalendarAgeDensityIndividualSample( 9, PP_fit_output, show_hpd_ranges = TRUE, show_unmodelled_density = TRUE) ## ----changepoint_locations_AD, out.width="100%"------------------------------- PlotPosteriorChangePoints(PP_fit_output, plot_cal_age_scale = "AD") ## ----find_meanrate, out.width="100%"------------------------------------------ # Calculating 2 sigma (95.4%) intervals on posterior mean occurrence rate PP_posterior_mean_rate_2sigma <- FindPosteriorMeanRate( PP_fit_output, calendar_age_sequence = seq(300, 500, by = 1), interval_width = "2sigma") # Look at posterior mean with 2 sigma probability interval head(PP_posterior_mean_rate_2sigma) ## ----plot_conditionalmeanrate, out.width="100%"------------------------------- # Conditional on TWO internal changes in the occurrence rate, # Calculate and plot the posterior mean rate over time # (with its 2 sigma intervals) conditional_2_changes_posterior_mean_rate <- PlotPosteriorMeanRate( PP_fit_output, n_changes = 2) # here n_changes must have length one (i.e., a single number) # Look at conditional posterior mean with 2 sigma probability interval head(conditional_2_changes_posterior_mean_rate) ## ----plot_ind_realisation, out.width="100%"----------------------------------- # Choose some nice plotting colours (from Okabe-Ito) realisation_colours <- c("#000000", "#E69F00", "#56B4E9", "#009E73", "#F0E442" ) # Plot 5 random realisations from posterior PlotRateIndividualRealisation( PP_fit_output, n_realisations = 5, plot_realisations_colour = realisation_colours)