Bayesian Estimation of Univariate Probability Distributions
Description
For a given dataset, this function serves to approximate (using a Metropolis algorithm) the posterior distribution of the parameters for some specified parametric probability distribution.
Usage
1 2 3 
Arguments
model 
A 
iter 
The requested number of iterations  the Markov Chain's length. 
tun 
A tuning constant; value by which the covariance matrix of the multivariate normal proposal will be multiplied  see References. 
trans.list 
A 
start 
A vector of starting values for the algorithm. If

variance 
Covariance matrix of the multivariate normal proposal
distribution. If 
prior 
A function that corresponds to the joint prior distribution (see Example). Note that the prior distribution will be evaluated on the transformed parameter space(s). 
burn 
Burnin period (see References). 
uniroot.interval 
Default is c(100,100). This interval is used
by 
pass.down.to.C 
If 
Details
This function uses a single block Metropolis algorithm with
multivariate normal proposal. For this function to work properly, all
parameters should be defined on the real line  parameter
transformation(s) might be required. If trans.list
is not
specified, the function will assume that the parameter distributions
are all defined on the real line (i.e., function(x) x
will be
used for each parameter). If no prior distribution is provided, an
improper prior distribution  uniform on the interval )Inf,+Inf( 
will be used for all parameters (i.e., prior distribution proportional
to 1  function(x) 1).
In order to minimize the number of arguments for metropolis
,
the function automatically computes the inverse of trans.list
:
this suppresses the need for the user to provide both the
"inverse transformation" and the "transformation". However, problems
may occur, and it is why the user is allowed to alter
uniroot.interval
. Depending on the number of errors reported,
future versions of this package may end up requesting that a list for
both the "inverse transformation" and the "transformation" be provided
by the user.
A nice list of references is provided below for more information on topics such as: MCMC algorithms, tuning of MetropolisHastings algorithms, MCMC convergence diagnostics, the Bayesian paradigm ...
Value
rate 
MCMC acceptance rate. This value is computed before
applying the burnin; i.e., it is computed for 
total.time 
Total computation time. 
sims.all 
Array containing all iterations. 
sims 
Array containing iterations after burnin. 
input 
Inputted 
iter 
Number of iterations. 
prior 
Prior distribution. 
burn 
Integer corresponding to the number of iterations to be discarded  burnin period. 
M 
Parameter vector whose elements correspond to the parameter
values (on the scales specified by 
V 
Covariance matrix computed using, after removing the burnin
period, the joint posterior distribution of the parameters (on the
scales specified by 
References
Gelman, A., Carlin, J.B., Stern, H.S., and Rubin, D.B. (2004). Bayesian data analysis, 2nd edition, Chapman & Hall/CRC.
Carlin, B.P, and Louis, T.A. (2009). Bayesian methods for data analysis. Chapman & Hall/CRC.
Gamerman, D., and Lopes H.F. (2006). Markov Chain Monte Carlo: Stochastic simulation for Bayesian inference. 2nd edition, Chapman & Hall/CRC.
Gilks, W.R., Richardson, S., and Spiegelhalter, D.J. (1996). Markov Chain Monte Carlo in Practice. Chapman & Hall.
See Also
plot.metropolis
, mle
Examples
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31  ### These examples should be rerun with, e.g., iter > 2000.
data(yarns)
x < yarns$x
fit.x < mle(x,'gamma',c(.1,.1))
bayes.x.no.prior < metropolis(model=fit.x,iter=150,
trans.list=list(function(x) x,function(x) exp(x)))
plot(bayes.x.no.prior)
# examples of prior distributions (note that these prior distribution
# are specified for the transformated parameters;
# i.e., in this case, 'meanlog' > 'meanlog' and 'sdlog' > 'ln.sdlog')
# for the scale parameter only
prior.1 < function(x) dnorm(x[2],.8,.1)
# for both parameters (joint but independent in this case)
prior.2 < function(x) dunif(x[1],3.4,3.6)*dnorm(x[2],1,1)
bayes.x.prior.2 < metropolis(model=fit.x,iter=150,
trans.list=list(function(x) x,function(x) exp(x)),prior=prior.2)
plot(bayes.x.prior.2)
# Example where 'model' is not from the class 'mle'; i.e.
# both 'start' and 'variance' need to be specified!
#x < rweibull(5,2,1)
x < c(0.9303492,1.0894917,0.9628029,0.6145032,0.4756699)
# Here 'fit.x < mle(x,'weibull',c(.1,.1))' is not used,
model.x < list(x=x,dist='weibull')
# and an informative prior distribution is considered to ensure a proper posterior distribution
prior.x < function(x) dnorm(x[1],log(2),.1)*dnorm(x[2],log(1),.1)
trans.list.x < list(function(x) exp(x), function(x) exp(x))
bayes.x < metropolis(model=model.x,iter=150,prior=prior.x,trans.list=trans.list.x,
pass.down.to.C=TRUE,start=c(0,0),variance=diag(.1,2,2))
