posterior: MCMC Sampling of Posterior Distributions

Description

Constructing MCMC samples of prior and posterior distibutions for the location (with optional trend), scale and shape parameters under the gev, order statistics or point process models.

Usage

posterior(n, init, prior, lh = c("none", "gev", "gpd", "pp","os"), ..., psd,
    burn = 0, thin = 1)

Value

A matrix with 1+floor(n/thin)-burn rows. Row labels give the iteration numbers. Column labels give parameter names.

An attribute ar is also returned. This is a matrix containing acceptence rates in the first row (the number of proposals accepted divided by the number of iterations) and ``external rates'' in the second (the number of proposals that resulted in a zero likelihood, divided by the number of iterations).

Arguments

n: The run-length; the number of sampled vectors (excluding init).
init: Numeric vector of length three/four, giving the initial values for the chain, taken to be iteration zero.
prior: An object of class "evprior", constructed using prior.prob, prior.quant or prior.norm.
lh: A character string specifying the likelihood; either "gev" for gev, "gpd" for gpd, "os" for order statistics, "pp" for Poisson process or "none" for none (the default). The latter can be used to sample from the prior distribution.
...: Arguments to the likelihood. Should include data unless lh is "none". Should also include thresh and noy if lh is "pp". Should include the vector trend if a linear trend on the location is implemented. See pplik for details.
psd: A vector of length three/four containing standard deviations for proposal distributions.
burn: The burn-in period (an integer); the first burn iterations (including init) are excluded from the chain.
thin: The thinning interval (an integer); iteration \(k\) is stored only if \(k\) mod thin is zero (and if \(k\) greater than or equal to burn).

Details

See the user's guide.

Examples

Run this code

data(rainfall)
prrain <- prior.quant(shape = c(38.9, 7.1, 47), scale = c(1.5, 6.3, 2.6))

n <- 100 ; t0 <- c(50.8, 1.18, 0.65) ; s <- c(25, .35, .07) ; b <- 20
rn.prior <- posterior(n, t0, prrain, "none", psd = s, burn = b)

t0 <- c(43.2, 7.64, 0.32) ; s <- c(2, .2, .07)
rn.post <- posterior(n, t0, prrain, "pp", data = rainfall, thresh = 40,
  noy = 54, psd = s, burn = b)