GetEstimate_BAY {HydroPortailStats}

Bayesian estimation of a distribution

Description

Returns MCMC samples from the posterior distribution.

Usage

GetEstimate_BAY(
  y,
  dist,
  prior,
  par0,
  mult = 0.1,
  eps = 0.1,
  batch.length = 100,
  batch.n = 100,
  moverate.min = 0.1,
  moverate.max = 0.5,
  mult.down = 0.9,
  mult.up = 1.1
)

Arguments

y

numeric vector, data

dist

character, distribution name

prior

list of lists, prior distributions. For each parameter to be estimated, the prior is a list of the form pr=list(dist=..., par=...). See example below.

par0

numeric vector, initial parameter guess. You may use GetEstimate_ROUGH().

mult

numeric, initial jump standard deviations are set to mult * abs(par0)

eps

numeric, where par0 is zero, initial jump standard deviations are set to eps (to avoid jumps of size zero)

batch.length

integer, MCMC parameter: length of each non-adaptive batch

batch.n

integer, MCMC parameter: number of batches (= adaptation period). Total number of simulations is nsim=batch.n*batch.length

moverate.min

numeric in (0;1), MCMC parameter: lower bound for the desired move rate interval

moverate.max

numeric in (0;1), MCMC parameter: upper bound for the desired move rate interval

mult.down

numeric in (0;1), MCMC parameter: multiplication factor used to decrease jump size when move rate is too low.

mult.up

numeric (>1, avoid 1/mult.down), MCMC parameter: multiplication factor used to increase jump size when move rate is too high.

Value

A list with the following components:

x

numeric matrix nsim*length(x0), MCMC simulations

fx

numeric vector, corresponding values f(x)

Examples

y=c(9.2,9.5,11.4,9.5,9.4,9.6,10.5,11.1,10.5,10.4)
prior1=list(dist='FlatPrior',par=NULL)
prior2=list(dist='LogNormal',par=c(1,1))
prior3=list(dist='Normal',par=c(0,0.25))
prior=list(prior1,prior2,prior3)
par0=GetEstimate_ROUGH(y,'GEV')$par
mcmc=GetEstimate_BAY(y,'GEV',prior,par0,batch.length=50,batch.n=50)
graphicalpar=par(mfrow=c(2,3))
plot(mcmc$x[,1],type='l'); plot(mcmc$x[,2],type='l'); plot(mcmc$x[,3],type='l')
hist(mcmc$x[,1]); hist(mcmc$x[,2]); hist(mcmc$x[,3])
par(graphicalpar)