blr: Bayesian Linear Regression

Description

This function fits a Bayesian linear regression model using scale invariant prior.

Usage

blr(y, w, mcmc = list(), prior = list(), marginal.likelihood = TRUE)

Arguments

a vector of response values.

a vector or matrix giving covariates of dimension n times ndimw

mcmc

a list giving the MCMC parameters. The list includes the following integers (with default values in parentheses): nblow (1000) giving the number of MCMC in transition period, nskip (1) giving the thinning interval, smcmc (1000) giving the number of MCMC for analysis.

prior

a list giving the prior information. The list includes the following parameters (default values specify the non-informative prior): beta_m0 and beta_v0 giving the hyperparameters of the multivariate normal distribution for parametric part including intercept, sigma2_m0 and sigma2_v0 giving the prior mean and variance of the inverse gamma prior for the scale parameter of response.

marginal.likelihood

a logical variable indicating whether the log marginal likelihood is calculated.

Value

An object of class blm representing the Bayesian spectral analysis model fit. Generic functions such as print and fitted have methods to show the results of the fit.

The MCMC samples of the parameters in the model are stored in the list mcmc.draws and the posterior samples of the fitted values are stored in the list fit.draws. The output list also includes the following objects:

post.est

posterior estimates for all parameters in the model.

lmarg

log marginal likelihood.

rsquarey

correlation between $y$ and $\hat{y}$.

call

the matched call.

mcmctime

running time of Markov chain from system.time().

Details

This generic function fits a Bayesian linear regression model using scale invariant prior.

Let $y_i$ and $w_i$ be the response and the vector of parametric predictors, respectively. The model for regression function is as follows.

$$y_i = w_i^T\beta + \epsilon_i, ~ i=1,\ldots,n, $$ where the error terms $\{\epsilon_i\}$ are a random sample from a normal distribution, $N(0,\sigma^2)$.

The conjugate priors are assumed for $\beta$ and $\sigma$: $$\beta | \sigma \sim N(m_{0,\beta}, \sigma^2V_{0,\beta}), \quad \sigma^2 \sim IG\Big(\frac{r_{0,\sigma}}{2}, \frac{s_{0,\sigma}}{2}\Big)$$

Examples

Run this code

# NOT RUN {
	#####################
	# Simulated example #
	#####################

	# Simulate data
	  set.seed(1)

	  n <- 100
	  w <- runif(n)
	  y <- 3 + 2*w + rnorm(n, sd = 0.8)						

	  # Fit the model with default priors and mcmc parameters
	  fout <- blr(y = y, w = w)

	  # Summary
	  print(fout)

	  # Fitted values
	  fit <- fitted(fout) 

	  # Plots
	  plot(fout)
# }