mcmc: Makes a MCMC chain using the Metropolis-Hastings algorithm

In UCL/ADMUR: Ancient Demographic Modelling Using Radiocarbon

Makes a MCMC chain using the Metropolis-Hastings algorithm

Description

Generates a single Markov Chain Monte Carlo using the Metropolis-Hastings algorithm

Usage

mcmc(PDarray, 
						startPars, 
						type, 
						timeseries=NULL,
						N=30000, 
						burn=2000, 
						thin=5, 
						jumps=0.02)

Arguments

PDarray	A Probability Density Array of calibrated dates, generated by phaseCalibrator .
startPars	A vector of parameter values for the algorithm to start at. Suggest using the parameters found from a Maximum Likelihood search. Must have an odd length with values between 0 and 1 for an n-CPL model, or length 1 values between -0.1 and 0.1 for an exponential model.
type	Choose from the following currently available models. Composite models can be achieved using a vector of more than one type. For example, c('norm','power') will be a composite model, where the first two parameters are the mean and SD, the 3rd and 4th parameters determine the power distribution component, for example if modelling taphonomy.
timeseries	Only if 'type' = 'timeseries', a data frame should be provided containing columns x and y that define the timeseries as year and value respectively.
N	The total number of proposals made, before the chain is reduced for burn-in and thinning.
burn	The number of proposals made at the beginning of the chain to be disregarded as burn-in.
thin	Specifies the proportion of proposals to be disregarded (after burn-in), such that the default 5 will only keep every 5th proposal.
jumps	Vector that determines the size of the random jump between the last parameters and the proposed parameters. A smaller value gives a smaller jump. Different jump values can be provided for each parameter. This can be tuned by observing how well mixed each parameter is in the chain.

Details

Generates a single MCMC chain using the Metropolis-Hastings algorithm, printing to screen progress every 1000 proposals. An acceptance ratio of around 0.4 to 0.5 should be sought by adapting the arguments 'burn' and 'jumps'. If the acceptance ratio is too low try reducing jump. A larger dataset will typically require smaller jumps.

Value

Returns a list including:

res	A 2D matrix of free parameter values (between 0 and 1) in the chain, after burn-in and thinning.
all.pars	A 2D matrix of free parameter values (between 0 and 1) in the chain of all N proposals.
acceptance.ratio	The proportion of proposals (after burn-in) that are accepted.

Examples

# randomly sample calendar dates from the toy model
set.seed(12345) 
N <- 350
cal <- simulateCalendarDates(toy, N)

# Convert to 14C dates. 
age <- uncalibrateCalendarDates(cal, shcal20)
data <- data.frame(age = age, sd = 50, phase = 1:N, datingType = '14C')


# Calibrate each phase, taking care to restrict to the modelled date range
CalArray <- makeCalArray(shcal20, calrange = range(toy$year), inc = 5)
PD <- phaseCalibrator(data, CalArray, remove.external = TRUE)

# Run MCMC for a 3-CPL model (5 parameters)
chain <- mcmc(PDarray = PD, 
	startPars = rep(0.5,5), 
	type='CPL', 
	jumps = 0.02)	

# Run MCMC for a 3-CPL model with taphonomy (5 + 2 parameters)
chain <- mcmc(PDarray = PD, 
	startPars = c(rep(0.5,5),10000,-1.5), 
	type=c('CPL', 'power'), 
	jumps = 0.02)	

UCL/ADMUR documentation built on Sept. 14, 2023, 11:41 a.m.