gibbs.lm: Linear model Gibbs sampling

In tkmckenzie/ikde: Iterative Kernel Density Estimation

Description

Fits a linear model using Gibbs sampling and estimates marginal likelihood as in Chib (1995)

Usage

gibbs.lm(X, y, priors = list(), burn.iter = 1000, sample.iter = 1000)

Arguments

X	Matrix of input variables
y	Vector of output variables
priors	A named list of parameter priors; should include beta.prior.mean (vector), beta.prior.var (matrix), sigma.sq.prior.shape (scalar), and sigma.sq.prior.rate (scalar)
burn.iter	Number of warmup iterations
sample.iter	Number of sampling iterations

Details

Uses a standard formulation of a linear model from which a Gibbs sampler can be derived. Specifically, for a model of the form

β\sim N(μ_β, Σ_β)

σ^2\sim Γ(s_σ, r_σ)

y = Xβ + \varepsilon

\varepsilon\sim N≤ft(0, \frac{1}{√{σ^2}} I\right),

Gibbs sampling can be performed using the conditional distributions

β|σ^2, X, y\sim N(\tildeμ_β, \tildeΣ_β)

σ^2|β, X, y\sim Γ^{-1}≤ft(\frac{N}2 + s_σ, \frac{e'e}2 + r_σ\right),

where N is the number of observations and

\tildeΣ_β = \frac{X'X}{σ^2} + Σ_β^{-1}

\tildeμ_β = \tildeΣ_β ≤ft(\frac{X'y}{σ^2} + Σ_β^{-1}μ_β\right)

e = y - Xβ.

Value

Returns an list with the following elements

samples	Named list of samples from the posterior, with elements "beta" and "sigma.sq"
log.marginal	Estimate of the model's log-marginal-likelihood
priors	List of priors used for the model

Examples

data(lm.generated)

X <- lm.generated$X
y <- lm.generated$y

gibbs.fit <- gibbs.lm(X, y,
                      priors = list(beta.prior.mean = rep(0, 4),
                                    beta.prior.var = 100 * diag(4),
                                    sigma.sq.prior.shape = 1,
                                    sigma.sq.prior.rate = 1))

print(apply(gibbs.fit$samples$beta, 2, mean)) # [1] 3.181184 1.643960 4.480879 1.213804
print(mean(gibbs.fit$samples$sigma.sq)) # [1] 97.52314
print(gibbs.fit$log.marginal) # [1] -389.001

tkmckenzie/ikde documentation built on May 13, 2019, 9:53 p.m.