as.parameters: Extract parameters from mcmcComposite object

In HelpersMG: Tools for Environmental Analyses, Ecotoxicology and Various R Functions

Extract parameters from mcmcComposite object

Description

Take a mcmcComposite object and create a vector object with parameter value at specified iteration.
If index="best", the function will return the parameters for the highest likelihood. It also indicates at which iteration the maximum lihelihood has been observed.
If index="last", the function will return the parameters for the last likelihood.
If index="median", the function will return the median value of the parameter.
if index="quantile", the function will return the probs defined by quantiles parameter.
If index="mode", the function will return the mode value of the parameter based on Asselin de Beauville (1978) method.
index can also be a numeric value. It uses all the chains being concatanated.
This function uses the complete iterations available if total is TRUE. Is adaptation part is never used.

Usage

as.parameters(
  x = stop("A result obtained after a MCMC analysis must be given."),
  total = TRUE,
  index = "best",
  chain = "all",
  probs = c(0.025, 0.5, 0.975),
  silent = FALSE
)

Arguments

x	A mcmcComposite obtained as a result of MHalgoGen() function
total	If TRUE, does not use the thinned results.
index	At which iteration the parameters must be taken, see description.
chain	The chain in which to get parameters; "all" is for all chains.
probs	Quantiles to be returned, see description.
silent	If TRUE, does not print any information.

Value

A vector with parameters at maximum likelihood or index position

Author(s)

Marc Girondot [e-mail]

References

Asselin de Beauville J.-P. (1978). Estimation non paramétrique de la densité et du mode, exemple de la distribution Gamma. Revue de Statistique Appliquée, 26(3):47-70.

See Also

Other mcmcComposite functions: MHalgoGen(), as.mcmc.mcmcComposite(), as.quantiles(), merge.mcmcComposite(), plot.PriorsmcmcComposite(), plot.mcmcComposite(), setPriors(), summary.mcmcComposite()

Examples

## Not run: 
library(HelpersMG)
require(coda)
x <- rnorm(30, 10, 2)
dnormx <- function(data, x) {
 data <- unlist(data)
 return(-sum(dnorm(data, mean=x['mean'], sd=x['sd'], log=TRUE)))
}
parameters_mcmc <- data.frame(Density=c('dnorm', 'dlnorm'), 
                              Prior1=c(10, 0.5), 
                              Prior2=c(2, 0.5), 
                              SDProp=c(1, 1), 
                              Min=c(-3, 0), 
                              Max=c(100, 10), 
                              Init=c(10, 2), 
                              stringsAsFactors = FALSE, 
                              row.names=c('mean', 'sd'))
mcmc_run <- MHalgoGen(n.iter=1000, parameters=parameters_mcmc, data=x, 
                      likelihood=dnormx, n.chains=1, n.adapt=100, 
                      thin=1, trace=1)
plot(mcmc_run, xlim=c(0, 20))
plot(mcmc_run, xlim=c(0, 10), parameters="sd")
mcmcforcoda <- as.mcmc(mcmc_run)
#' heidel.diag(mcmcforcoda)
raftery.diag(mcmcforcoda)
autocorr.diag(mcmcforcoda)
acf(mcmcforcoda[[1]][,"mean"], lag.max=20, bty="n", las=1)
acf(mcmcforcoda[[1]][,"sd"], lag.max=20, bty="n", las=1)
batchSE(mcmcforcoda, batchSize=100)

# The batch standard error procedure is usually thought to 
# be not as accurate as the time series methods used in summary
summary(mcmcforcoda)$statistics[,"Time-series SE"]
summary(mcmc_run)
as.parameters(mcmc_run)
lastp <- as.parameters(mcmc_run, index="last")
parameters_mcmc[,"Init"] <- lastp

# The n.adapt set to 1 is used to not record the first set of parameters
# then it is not duplicated (as it is also the last one for 
# the object mcmc_run)
mcmc_run2 <- MHalgoGen(n.iter=1000, parameters=parameters_mcmc, data=x, 
                       likelihood=dnormx, n.chains=1, n.adapt=1, 
                       thin=1, trace=1)
mcmc_run3 <- merge(mcmc_run, mcmc_run2)

####### no adaptation, n.adapt must be 0
parameters_mcmc[,"Init"] <- c(mean(x), sd(x))
mcmc_run3 <- MHalgoGen(n.iter=1000, parameters=parameters_mcmc, data=x, 
                       likelihood=dnormx, n.chains=1, n.adapt=0, 
                       thin=1, trace=1)
# With index being median, it returns the median value for each parameter
as.parameters(mcmc_run3, index="median")
as.parameters(mcmc_run3, index="mode")
as.parameters(mcmc_run3, index="best")
as.parameters(mcmc_run3, index="quantile", probs=0.025)
as.parameters(mcmc_run3, index="quantile", probs=0.975)
as.parameters(mcmc_run3, index="quantile", probs=c(0.025, 0.975))

## End(Not run)

HelpersMG documentation built on April 4, 2025, 4:40 a.m.