normGibbs {Bolstad2} | R Documentation |
Draw a sample from a posterior distribution of data with an unknown mean and variance using Gibbs sampling
Description
normGibbs draws a Gibbs sample from the posterior distribution of the
parameters given the data fron normal distribution with unknown mean and
variance. The prior for given
is prior mean
and prior variance
. That means
is the 'equivalent
sample size.' The prior distribution of the variance is
times an
inverse chi-squared with
degrees of freedom. The joint prior is
the product
.
Usage
normGibbs(y, steps = 1000, type = "ind", ...)
Arguments
y |
A vector containing the data |
steps |
The number of iterations of Gibbs sampling to carry out |
type |
Either 'ind' for sampling from an independent conjugate prior or 'joint' for sampling from a joint conjugate prior. 'i' and 'j' can be used as compact notation |
... |
If type = 'ind' then the user can specify the prior for
|
Value
A data frame containing three variables
[1,] | mu | a sample from the posterior distribution of the mean |
[2,] | sig | a sample from the posterior distribution of the standard deviation |
[3,] | mu | a sample from the posterior distribution of the variance = sig^2 |
Author(s)
James M. Curran
Examples
## firstly generate some random data
mu = rnorm(1)
sigma = rgamma(1,5,1)
y = rnorm(100, mu, sigma)
## A \eqn{N(10,3^2)} prior for \eqn{\mu} and a 25 times inverse chi-squared
## with one degree of freedom prior for \eqn{\sigma^2}
MCMCSampleInd = normGibbs(y, steps = 5000, priorMu = c(10,3),
priorVar = c(25,1))
## We can also use a joint conjugate prior for \eqn{\mu} and \eqn{\sigma^2}.
## This will be a \emph{normal}\eqn{(m,\sigma^2/n_0)} prior for \eqn{\mu} given
## the variance \eqn{\sigma^2}, and an \eqn{s0} times an \emph{inverse
## chi-squared} prior for \eqn{\sigma^2}.
MCMCSampleJoint = normGibbs(y, steps = 5000, type = 'joint',
priorMu = c(10,3), priorVar = c(25,1))
## Now plot the results
oldPar = par(mfrow=c(2,2))
plot(density(MCMCSampleInd$mu),xlab=expression(mu), main =
'Independent')
abline(v=mu)
plot(density(MCMCSampleInd$sig),xlab=expression(sig), main =
'Independent')
abline(v=sigma)
plot(density(MCMCSampleJoint$mu),xlab=expression(mu), main =
'Joint')
abline(v=mu)
plot(density(MCMCSampleJoint$sig),xlab=expression(sig), main =
'Joint')
abline(v=sigma)