phenology (version 4.0.4)

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

result

An object obtained after a SearchR fit

n.iter

Number of iterations for each step

parametersMCMC

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

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))