probit_regression: Gibbs Sampler for Probit Regression

In dcbdan/s525: Implementation of common Bayesian models

Description

The probit regression model is given by

y_i|π_i \sim Bernoulli(π_i),

π_i = Φ(x_i^tβ),

β_j \sim N(0, A^2).

where Φ is given by the gaussian CDF.

The implemented parameter-expanded model is given by

y_i^* = x_i^tβ + ε_i,

ε_i \sim N(0, 1),

y_i = I(y_i^* > 0).

The full conditional distributions are given by

[β|y, y^*, X] \sim N(Q^{-1}l, Q^{-1})

where Q = X^tX + A^{-2}I and l = X^ty^*,

[y_i^*|y_i, β, X] \sim sgn(y_i, y_i^*)N(x_i^tβ, 1)

where sqn is 1 if both arguments are of the same sign and zero otherwise.

If A is a vector, Q = X^tX + diag(A). A flat prior on the intercept at position 1 is thus given by A = c(Inf, rep(3, p)).

Usage

probit_regression <- function(
  Y,
  X,
  niter,
  A = 3,
  init = NULL)

Arguments

Y	n by 1 vector of ones and zeros
X	n by p predictor matrix
niter	number of gibbs sampling iterations
A	A parameter. The default value is chosen to provide a reasonable range of π. If A is a vector, then different variances are given for each intercept.
init	Initial starting values for beta. If NULL, beta is set to zero.

Details

This function returns a niter x p matrix of values where p is the second dimension of the predictor matrix X. The returned matrix contains all generated values of the gibbs sampling markov chain.

Examples

  library(LearnBayes)
  library(s525)

  data(donner)
  y = donner$survival
  X = cbind(1, donner$age, donner$male)
  niter <- 2000

  gibbs_results <- probit_regression(y, X, niter)

dcbdan/s525 documentation built on May 19, 2019, 10:48 p.m.