phenology (version 3.72)

phenology_MHmcmc: Run the Metropolis-Hastings algorithm for data

Description

Run the Metropolis-Hastings algorithm for data. Deeply modified from a MCMC script by Olivier Martin (INRA, Paris-Grignon). The number of iterations is n.iter+n.adapt+1 because the initial likelihood is also displayed. I recommend that thin=1 because the method to estimate SE uses resampling. As initial point is maximum likelihood, n.adapt = 0 seems a good solution.

Usage

phenology_MHmcmc(result = stop("An output from fit_phenology() must be provided"),
  n.iter = 10000,
  parametersMCMC = stop("A model generated with phenology_MHmcmc_p() must be provided"),
  n.chains = 4, n.adapt = 0, thin = 1, trace = FALSE)

Arguments

n.iter

Number of iterations for each step

parametersMCMC

A set of parameters used as initial point for searching with information on priors

result

An object obtained after a SearchR fit

n.chains

Number of replicates

n.adapt

Number of iterations before to store outputs

thin

Number of iterations between each stored output

trace

True or False, shows progress

Value

• A list with resultMCMC being mcmc.list object, resultLnL being likelihoods and parametersMCMC being the parameters used

Details

phenology_MHmcmc runs the Metropolis-Hastings algorithm for data (Bayesian MCMC)

Examples

library(phenology)
data(Gratiot)
# Generate a formatted list named data_Gratiot
data_Gratiot<-add_phenology(Gratiot, name="Complete",
    reference=as.Date("2001-01-01"), format="%d/%m/%Y")
# Generate initial points for the optimisation
parg<-par_init(data_Gratiot, parametersfixed=NULL)
# Run the optimisation
result_Gratiot<-fit_phenology(data=data_Gratiot,
		parametersfit=parg, parametersfixed=NULL, trace=1)
# Generate set of priors for Bayesian analysis
pmcmc <- phenology_MHmcmc_p(result_Gratiot, accept = TRUE)
result_Gratiot_mcmc <- phenology_MHmcmc(result = result_Gratiot, n.iter = 10000,
parametersMCMC = pmcmc, n.chains = 1, n.adapt = 0, thin = 1, trace = FALSE)
# Get standard error of parameters
summary(result_Gratiot_mcmc)
# Make diagnostics of the mcmc results using coda package
mcmc <- as.mcmc(result_Gratiot_mcmc)
require(coda)
heidel.diag(mcmc)
raftery.diag(mcmc)
autocorr.diag(mcmc)
acf(mcmc[[1]][,"LengthB"], lag.max=200, bty="n", las=1)
acf(mcmc[[1]][,"Max_Gratiot"], lag.max=50, bty="n", las=1)
batchSE(mcmc, batchSize=100)
# The batch standard error procedure is usually thought to
# be not as accurate as the time series methods used in summary
summary(mcmc)$statistics[,"Time-series SE"]
plot(result_Gratiot_mcmc, parameters=3, las=1, xlim=c(-10, 300))