R: Draw a sample from a posterior distribution of data with an...
normGibbs
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 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)
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/normGibbs.Rd_%03d_medium.png", width=480, height=480)
> ### Name: normGibbs
> ### Title: Draw a sample from a posterior distribution of data with an
> ### unknown mean and variance using Gibbs sampling
> ### Aliases: normGibbs
>
> ### ** 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))
N mean stdev sterr min q1 med
mu 5000 -1.238207 0.8580427 0.012134557 -5.405055 -1.824605 -1.237459
sig 5000 8.501562 0.6136423 0.008678212 6.661707 8.079131 8.474729
var 5000 72.653034 10.6291424 0.150318773 44.378339 65.272352 71.821035
q3 max
mu -0.6586607 2.352413
sig 8.8839139 13.939501
var 78.9239260 194.309674
>
>
> ## 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))
N mean stdev sterr min q1 med
mu 5000 -1.232931 0.8328797 0.011778698 -4.119827 -1.783157 -1.246462
sig 5000 8.480864 0.6166467 0.008720701 6.663897 8.053846 8.434067
var 5000 72.305228 10.6338238 0.150384978 44.407519 64.864434 71.133485
q3 max
mu -0.6662507 1.803212
sig 8.8673604 11.708424
var 78.6300803 137.087198
>
> ## 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)
>
>
>
>
>
>
> dev.off()
null device
1
>
providing 
.
