normGibbs: Draw a sample from a posterior distribution of data with an...

View source: R/normGibbs.R

normGibbsR 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 var is prior mean m0 and prior variance var/n0 . That means n0 is the 'equivalent sample size.' The prior distribution of the variance is s0 times an inverse chi-squared with kappa0 degrees of freedom. The joint prior is the product g(var)g(mu|var).

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 μ with a parameter priorMu which can either be a single number m0, or m0 and n0. if m0 and n0 are not specified then m0 and n0 are 0 by default. The user can also specify priorVar, which if given, must be a vector with two elements s0 and kappa0. If s0 and kappa0 are not given then they are zero by default. If type = 'joint' then priorMu must be a vector of length two with elements m0 and sd0. The user can also specify priorVar, which if given, must be a vector with two elements s0 and kappa0. If s0 and kappa0 are not given then they are zero by default.

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)



Bolstad2 documentation built on April 11, 2022, 5:08 p.m.

Related to normGibbs in Bolstad2...