findStartingValues: Find Starting Values for the Metropolis Algorithm
In ClementCalenge/metroponcfs: Metropolis algorithm for model fitting

Description Usage Arguments Details Value Author(s) See Also Examples

Find Starting Values for the Metropolis Algorithm

1
2
3

findStartingValues(support, logposterior, lidat, multiple = list(),
  nrand = 10, method = c("onebest", "dispersed", "all"), ndispersed = 3,
  info = TRUE)

`support`	named list with one element per parameter (the names of the list correspond to the name of the parameters). Each element should be a vector of length 2 indicating the limits of the support of the parameters (i.e. the min and max values in which the starting values will be searched).
`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 (see the help page of `GeneralSingleMetropolis`). Note that the argument `ctrl` is not used by the function here.
`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
`multiple`	a named list describing which parameters are vectors of parameters. For example, if a model relies on a single parameter `"a"` and two vectors of parameters `"theta1"` and `"theta2"` of length 5 and 6 respectively, the argument `multiple` should take the value: `list(theta1=5, theta2=6)`. Note that single parameters are not indicated in this list.
`nrand`	the number of starting values to generate before a result is returned (see details).
`method`	the method to be used by the function (see details).
`ndispersed`	integer. When `method="dispersed"`, indicates how many dispersed starting values should be returned. Should be greater than `nrand`.
`info`	whether information on the generation process should be displayed to the user.

findStartingValues allows to find a list of starting values for the Metropolis algorithm for the MCMC fitting of Bayesian models.

If method=="onebest", a named list of starting values for the parameters, that can be used directly as an argument to GeneralSingleMetropolis. For other methods, a list of lists of starting values.

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

GeneralSingleMetropolis

## We generate data
x1 <- rnorm(100)
x2 <- rnorm(100)
x3 <- rnorm(100)
y <- 1+x1+0.5*x2+1.5*x3+rnorm(100, sd=0.5)

## Our data
lid <- list(x=cbind(x1,x2,x3), y=y)

## 5 parameters: intercept + 3 slopes + residual variance
## For the sake of demonstration, we put the three slopes in
## a vector theta for the estimation.
## The list of starting values should therefore contain: (i) theta,
## (ii) intercept, (iii) residprec
##
## The function for the posterior is:
lpost <- function(par, lidat, ctrl)
{
     ## prior:
     parb <- par
     parb$residprec <- NULL
     lprior <- sum(dnorm(unlist(parb), mean=0, sd=100)) ## vague prior
     ## (uniform prior for the residual precision, we do not add this constant)

     mu <- lidat$x%*%par$theta+par$intercept
     ## log-posterior
     return(sum(lprior + dnorm(lidat$y, mean=mu, sd = 1/sqrt(par$residprec), log=TRUE)))
}

##  We define large support for all parameters
support <- list(theta=c(-10,10), intercept=c(-10,10), residprec=c(0.0001, 10))
## only theta is multiple (3 elements), so:
multiple <- list(theta=3)

## Default method: finds the best starting value among 1000 sampled values
## The plot shows the distribution of the sampled posterior (up to a constant)
## the best value is returned in sv
(sv <- findStartingValues(support, lpost, lid, multiple, nrand=1000))

## Same, but keeps 4 most dispersed starting values among 100 generated values:
(sv2 <- findStartingValues(support, lpost, lid, multiple, nrand=100,
                           method="dispersed", ndispersed=4))