GeneralSingleMetropolis: Implements one Markov chain with the Metropolis algorithm
In ClementCalenge/metroponcfs: Metropolis algorithm for model fitting

Description Usage Arguments Details Value Note Author(s) Examples

Implements one Markov chain with the Metropolis algorithm

GeneralSingleMetropolis(parInit, lidat, logposterior, lipum,
  listUpdating = NULL, liopt = NULL, nrepet = 999, nburnin = 10,
  thinPar = 1, thinAcc = thinPar, composition = c("fixed", "random"),
  verbose = TRUE, saveResults = TRUE, saveEvery = 1000,
  fileSave = "saveMetrop.Rdata", debug = FALSE,
  optionsDebug = list(showSeed = FALSE, showPar = TRUE))

restartGSM(resuParms, another = 1000)

defaultListUGSM(parInit)

## S3 method for class 'GSMetrop'
print(x, ...)

reloadGSM(fileSave)

`parInit`	named list with one element per parameter (the names of the list correspond to the name of the parameters), containing the initial values for the parameters
`lidat`	a named list containing the data (the names of the list correspond to the name of the variables used in logposterior, combined with the parameters), required for the calculation of the posterior.
`logposterior`	a function for the calculation of the log posterior for the current model. This function must rely only on the data provided in `lidat`. This function must have three parameters: (i) `par` is the list of parameters, (ii) `lidat` is the list containing the data, and (iii) `ctrl` is a named list with two elements: (i) an element `namePar` is a character string containing the name of the parameter in `par` that is updated, (ii) if this parameter is a vector where each component is updated in turn (e.g. when `MultipleIndependentSteps` is the updating mechanism), `iter` indicates the element of this vector that is updated (this element is ignored in other cases). It is possible to define a log-posterior that does not make use of `ctrl` (i.e., which returns the same value whatever the parameter that is updated), but even in this case, the function should have an argument named `ctrl`. Note that the argument `ctrl` is essentially useful when the model contains many parameters (e.g. overdispersion residuals).
`lipum`	a named list containing the parameters of the updating mechanisms (see the help page of the updating mechanisms for a description of the required parameters). The names of the list correspond to the name of the parameters of the model (i.e. `names(parInit)`).
`listUpdating`	named list with one element parameter (the names of the list correspond to the name of the parameters), each one containing one of the following character strings: (i) "mis", which implements the function `MultipleIndependentSteps` as updating function for the corresponding parameter, (ii) "sis", which implements the function `singleStep` as updating function for the corresponding parameter, (iii) "mns", which implements the function `MultiNormalStep`. If `NULL`, the function `defaultListUGSM` is used internally to generate a default value ("mis" or "sis").
`liopt`	optionally, a vector containing the name of the parameters for which the log-posterior is optimized (i.e. the function logposterior returns a value that is proportional to the log-posterior, and not the log-poserior itself, see the help page of `singleStep`).
`nrepet`	an integer giving the number of iterations required for the MCMC algorithm.
`nburnin`	an integer giving the number of burnin iterations required
`thinPar`	an integer giving the number of steps separating the storage of two parameter vectors generated by the algorithm
`thinAcc`	an integer giving the number of steps separating the storage of the vector describing whether the proposal has been accepted or not
`composition`	character. Determines in which order the parameters are updated. Either `"fixed"` (the default), in which case the parameters are updated in the order provided in `parInit`, or `"random"`, in which case the parameters are updated in a random order at each step.
`verbose`	logical. Whether information should be displayed to the user.
`saveResults`	logical. Whether the results should be saved in a file.
`saveEvery`	integer. When should the program save the list of generated values (for backup)?
`fileSave`	character string. The name of the Rdata file used to save the results.
`debug`	logical. Whether the chain should be debugged (in this case, the whole state of the chain is stored in a dump file, storing a list with one element per iteration, storing the whole process). Debug mode is not possible for the function `restartGSM`.
`optionsDebug`	list with two elements named `showPar` and `showSeed` indicating what to show in the dump file at each iteration (at each iteration, the seed of the random number generator can be saved, as well as the value of the list of parameters, so that the user can focus on a particular step without bothering of the other steps).
`resuParms`	an object of class GSMetrop.
`another`	an integer giving the number of iterations required for the MCMC algorithm.
`x`	an object of class GSMetrop.
`...`	additional arguments to be passed to or from other functions

GeneralSingleMetropolis implements the Metropolis algorithm for the MCMC fitting of Bayesian models. restartGSM can be used to restart a model for another iterations. reloadGSM can be used to return the GSMetro object saved in a file by the backup mechanism of GeneralSingleMetropolis. defaultListUGSM can be used to generate default values for listUpdating (defaulting to "sis" for parameters of length 1 and to "mis" for vectors of parameters). This list can then be manually altered (e.g. changing the updating mechanism to "mns" for some parameters).

A list of class "GSMetrop", where each element corresponds to one parameter, and where each element is a matrix containing the generated values as rows.

The debug mode produces a temporary file with a name starting by dumpmetrop***.R, which contains the R code corresponding to the state of the chain at each step (old value of the parameter, proposed value, etc.). This file contains all the elements required to understand the state of the chain. At the top of the file, the function includes R code that allows to explore the chain in another session (i.e. the code in this file does not rely on the global environment). The arguments of the functions are restored in the environment at the beginning of the file.

Clement Calenge, clement.calenge@oncfs.gouv.fr

## Generates a dataset:
x <- rnorm(100)
y <- 1+ 2*x + rnorm(100, sd = 0.1)

## Linear regression: estimate a, b, sigma in
## y = a + b*x + epsilon
## epsilon ~ dnorm(0, sigma)
## We will model lsigma=log(sigma) to ensure it is >0

## Log-posterior
logp <- function(par, lidat, ctrl)
{
   ## Current values of the parameters
   a <- par$a
   b <- par$b
   sigma <- exp(par$lsigma) ## sigma=exp(log(sigma))

   ## log-Prior:
   lprior <- dnorm(a, mean=0, sd=100, log=TRUE)
   lprior <- lprior + dnorm(b, mean=0, sd=100, log=TRUE)
   ## uniform prior on log-sigma

   ## log-likelihood
   llike <- sum(dnorm(lidat$y, mean=lidat$x*b+a, sd=sigma, log=TRUE))

   ## log-posterior
   return(llike+lprior)
}

## Elements required for the fit
parinit <- list(a=0, b=1, lsigma=0) ## initial values

lipum <- list(a=0.1, b=0.1, lsigma=1) ## test standard deviation of the proposals

lidat <- list(x=x, y=y) ## the data

## Not run: 
gsm <- GeneralSingleMetropolis(parinit, lidat, logp, lipum,
                               nrepet=10000, saveResults = FALSE)
gsm

## Plot the chain
plot(gsm)

## We increase the burnin
gsm2 <- increaseBurnin(gsm, newBurnin=200)

## Ok
plot(gsm2)

## show a summary of the results
summary(gsm2, parameters=TRUE, accept=TRUE)

## Try to use a multinormal updating for a and b
## we change the logposterior, and define theta
## as c(a,b)
logp2 <- function(par, lidat, ctrl)
{
   ## Current values of the parameters
   theta <- par$theta
   a <- theta[1]
   b <- theta[2]
   sigma <- exp(par$lsigma)

   ## log-Prior:
   lprior <- dnorm(a, mean=0, sd=100, log=TRUE)
   lprior <- lprior + dnorm(b, mean=0, sd=100, log=TRUE)
   ## uniform prior on sigma

   ## log-likelihood
   llike <- sum(dnorm(lidat$y, mean=lidat$x*b+a, sd=sigma, log=TRUE))

   ## log-posterior
   return(llike+lprior)
}

## We get the last value of the parameters generated in gsm2
pr <- getParameterVector(gsm2, nrow(gsm2[[1]]))

## We change the variables
parinit2 <- list(theta=c(pr$a, pr$b), lsigma=pr$lsigma) ## initial values
listUpdating2 <- list(theta="mns", lsigma="sis") ## now, theta is updated with mns
lipum2 <- list(theta=cov(cbind(gsm2$a, gsm2$b)), lsigma=1)

## And we start again, saving the results in ficxample231.Rdata
gsm3 <- GeneralSingleMetropolis(parinit2, lidat, logp2, lipum2, listUpdating2,
                                nrepet=10000, fileSave="ficxample231.Rdata")

gsm3

plot(gsm3)

summary(gsm3, parameters=TRUE, accept=TRUE)

## Demonstrate how to continue the iterations:
gsm4 <- restartGSM(gsm3, another=10000)
gsm4

## demonstrates how to reload a file:
gsm5 <- reloadGSM("ficxample231.Rdata")

## housekeeping
file.remove("ficxample231.Rdata")

## End(Not run)