mcmc {ADMUR}

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, 
						taphonomy=FALSE, 
						taph.min=c(0,-3), 
						taph.max=c(20000,0), 
						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 'CPL', 'exp', 'norm', 'cauchy', 'sine', 'uniform', 'logistic', 'power'.

taphonomy

If TRUE, the last two parameters determine the taphonomic loss rate.

taph.min

Minimum values for the prior range of taphonomic parameters (b,c).

taph.max

Maximum values for the prior range of taphonomic parameters (b,c).

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='CPL', 
	taphonomy=TRUE, 
	jumps = 0.02)