View source: R/phenology_MHmcmc.R
phenology_MHmcmc | R Documentation |
Run the Metropolis-Hastings algorithm for data.
The number of iterations is n.iter+n.adapt+1 because the initial likelihood is also displayed.
I recommend thin=10.
If initial point is maximum likelihood, n.adapt = 0 is a good solution.
The parameters intermediate and filename are used to save intermediate results every 'intermediate' iterations (for example 1000). Results are saved in a file of name filename.
The parameter previous is used to indicate the list that has been save using the parameters intermediate and filename. It permits to continue a mcmc search.
These options are used to prevent the consequences of computer crash or if the run is very very long and computer processes are time limited.
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 = 1,
n.adapt = 1000,
thin = 1,
trace = FALSE,
traceML = FALSE,
adaptive = TRUE,
adaptive.lag = 500,
adaptive.fun = function(x) {
ifelse(x > 0.234, 1.3, 0.7)
},
intermediate = NULL,
filename = "intermediate.Rdata",
previous = NULL
)
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 or period, shows progress |
traceML |
TRUE or FALSE to show ML |
adaptive |
Should an adaptive process for SDProp be used |
adaptive.lag |
Lag to analyze the SDProp value in an adaptive content |
adaptive.fun |
Function used to change the SDProp |
intermediate |
Period for saving intermediate result, NULL for no save |
filename |
If intermediate is not NULL, save intermediate result in this file |
previous |
Previous result to be continued. Can be the filename in which intermediate results are saved. |
phenology_MHmcmc runs the Metropolis-Hastings algorithm for data (Bayesian MCMC)
A list with resultMCMC being mcmc.list object, resultLnL being likelihoods and parametersMCMC being the parameters used
Marc Girondot
Other Phenology model:
AutoFitPhenology()
,
BE_to_LBLE()
,
Gratiot
,
LBLE_to_BE()
,
LBLE_to_L()
,
L_to_LBLE()
,
MarineTurtles_2002
,
MinBMinE_to_Min()
,
adapt_parameters()
,
add_SE()
,
add_phenology()
,
extract_result()
,
fit_phenology()
,
likelihood_phenology()
,
logLik.phenology()
,
map_Gratiot
,
map_phenology()
,
par_init()
,
phenology()
,
phenology2fitRMU()
,
phenology_MHmcmc_p()
,
plot.phenology()
,
plot.phenologymap()
,
plot_delta()
,
plot_phi()
,
print.phenology()
,
print.phenologymap()
,
print.phenologyout()
,
remove_site()
,
result_Gratiot
,
result_Gratiot1
,
result_Gratiot2
,
result_Gratiot_Flat
,
result_Gratiot_mcmc
,
summary.phenology()
,
summary.phenologymap()
,
summary.phenologyout()
## Not run:
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, fixed.parameters=NULL)
# Run the optimisation
result_Gratiot <- fit_phenology(data=data_Gratiot,
fitted.parameters=parg, fixed.parameters=NULL)
# 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))
## End(Not run)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.