either a sample from the posterior density or the values over which the the posterior CDF is specified
cdf
the values of the CDF, F(θ) =
int_{-∞}^{θ}f(t).df where f(t) is the PDF. This
only needs to be specified if a numerically specified posterior is
being used
conf
the desired 'confidence' level
type
the type of interval to return, 'lower' = one sided lower
bound, 'two-sided' = two - sided, or 'upper' = one sided upper
bound. It is sufficient to use 'l','t' or 'u'
Details
This function uses linear interpolation to calculate bounds for points
that may not be specified by CDF
Value
a list containing the elements lower.bound, uppper.bound or both
depending on type
Examples
## commands for calculating a numerical posterior CDF.
## In this example, the likelihood is proportional to
## eqn{\theta^{3/2}\times exp(-\theta/4)} and a N(6, 9) prior is used.
theta <- seq(from = 0.001, to = 40, by = 0.001)
prior <- dnorm(theta,6,3)
ppnLike <- theta^1.5*exp(-theta/4)
ppnPost <- prior*ppnLike
scaleFactor <- sintegral(theta, ppnPost)$int
posterior <- ppnPost/scaleFactor
cdf <- sintegral(theta, posterior)$y
ci<-credInt(theta, cdf)
par(mfrow=c(2,2))
plot(prior ~ theta, type = 'l', main = "Prior N(6, 9)")
plot(ppnLike ~ theta, type = 'l', main = "Proportional likelihood")
plot(posterior ~ theta, type = 'l', main = "Posterior")
abline(v=c(unlist(ci)))
## Use an inverse method to take a random sample of size 1000
## from the posterior
suppressWarnings(Finv<-approxfun(cdf,theta))
thetaSample<-Finv(runif(1000))
ci<-credInt(thetaSample)
Results
R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> library(Bolstad2)
> png(filename="/home/ddbj/snapshot/RGM3/R_CC/result/Bolstad2/credInt.Rd_%03d_medium.png", width=480, height=480)
> ### Name: credInt
> ### Title: Calculate a credible interval from a numerically specified
> ### posterior CDF or from a sample from the posterior
> ### Aliases: credInt
>
> ### ** Examples
>
> ## commands for calculating a numerical posterior CDF.
> ## In this example, the likelihood is proportional to
> ## eqn{\theta^{3/2}\times exp(-\theta/4)} and a N(6, 9) prior is used.
> theta <- seq(from = 0.001, to = 40, by = 0.001)
> prior <- dnorm(theta,6,3)
> ppnLike <- theta^1.5*exp(-theta/4)
> ppnPost <- prior*ppnLike
> scaleFactor <- sintegral(theta, ppnPost)$int
> posterior <- ppnPost/scaleFactor
> cdf <- sintegral(theta, posterior)$y
> ci<-credInt(theta, cdf)
Credible interval is : (2.02114790458962,11.4246333371552)
> par(mfrow=c(2,2))
> plot(prior ~ theta, type = 'l', main = "Prior N(6, 9)")
> plot(ppnLike ~ theta, type = 'l', main = "Proportional likelihood")
> plot(posterior ~ theta, type = 'l', main = "Posterior")
> abline(v=c(unlist(ci)))
>
> ## Use an inverse method to take a random sample of size 1000
> ## from the posterior
> suppressWarnings(Finv<-approxfun(cdf,theta))
> thetaSample<-Finv(runif(1000))
> ci<-credInt(thetaSample)
Credible interval is (2.2300076598049,11.0198564940339)
>
>
>
>
>
> dev.off()
null device
1
>