Gibbs Sampler for Univariate Regression
runiregGibbs implements a Gibbs Sampler to draw from posterior of a univariate regression with a conditionally conjugate prior.
runiregGibbs(Data, Prior, Mcmc)
Data |
list(y, X) |
Prior |
list(betabar, A, nu, ssq) |
Mcmc |
list(sigmasq, R, keep, nprint) |
y = Xβ + e with e ~ N(0, σ^2)
β ~ N(betabar, A^{-1})
σ^2 ~ (nu*ssq)/χ^2_{nu}
Data = list(y, X)
y: |
n x 1 vector of observations |
X: |
n x k design matrix |
Prior = list(betabar, A, nu, ssq) [optional]
betabar: |
k x 1 prior mean (def: 0) |
A: |
k x k prior precision matrix (def: 0.01*I) |
nu: |
d.f. parameter for Inverted Chi-square prior (def: 3) |
ssq: |
scale parameter for Inverted Chi-square prior (def: var(y))
|
Mcmc = list(sigmasq, R, keep, nprint) [only R required]
sigmasq: |
value for σ^2 for first Gibbs sampler draw of β|σ^2 |
R: |
number of MCMC draws |
keep: |
MCMC thinning parameter -- keep every keepth draw (def: 1) |
nprint: |
print the estimated time remaining for every nprint'th draw (def: 100, set to 0 for no print)
|
A list containing:
betadraw |
R x k matrix of betadraws |
sigmasqdraw |
R x 1 vector of sigma-sq draws |
Peter Rossi, Anderson School, UCLA, perossichi@gmail.com.
For further discussion, see Chapter 3, Bayesian Statistics and Marketing by Rossi, Allenby, and McCulloch.
http://www.perossi.org/home/bsm-1
if(nchar(Sys.getenv("LONG_TEST")) != 0) {R=1000} else {R=10}
set.seed(66)
n = 200
X = cbind(rep(1,n), runif(n))
beta = c(1,2)
sigsq = 0.25
y = X%*%beta + rnorm(n,sd=sqrt(sigsq))
out = runiregGibbs(Data=list(y=y, X=X), Mcmc=list(R=R))
cat("Summary of beta and Sigmasq draws", fill=TRUE)
summary(out$betadraw, tvalues=beta)
summary(out$sigmasqdraw, tvalues=sigsq)
## plotting examples
if(0){plot(out$betadraw)}Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.