update_Sigma: Update error term covariance matrix of multiple linear...

In RprobitB: Bayesian Probit Choice Modeling

Update error term covariance matrix of multiple linear regression

Description

This function updates the error term covariance matrix of a multiple linear regression.

Usage

update_Sigma(kappa, E, N, S)

Arguments

kappa	The degrees of freedom (a natural number greater than J-1) of the Inverse Wishart prior for Sigma. Per default, kappa = J + 1.
E	The scale matrix of dimension J-1 x J-1 of the Inverse Wishart prior for Sigma. Per default, E = diag(J - 1).
N	The draw size.
S	A matrix, the sum over the outer products of the residuals (ε_n)_{n=1,…,N}.

Details

This function draws from the posterior distribution of the covariance matrix Σ in the linear utility equation

U_n = X_nβ + ε_n,

where U_n is the (latent, but here assumed to be known) utility vector of decider n = 1,…,N, X_n is the design matrix build from the choice characteristics faced by n, β is the coefficient vector, and ε_n is the error term assumed to be normally distributed with mean 0 and unknown covariance matrix Σ. A priori we assume the (conjugate) Inverse Wishart distribution

Σ \sim W(κ,E)

with κ degrees of freedom and scale matrix E. The posterior for Σ is the Inverted Wishart distribution with κ + N degrees of freedom and scale matrix E^{-1}+S, where S = ∑_{n=1}^{N} ε_n ε_n' is the sum over the outer products of the residuals (ε_n = U_n - X_nβ)_n.

Value

A matrix, a draw from the Inverse Wishart posterior distribution of the error term covariance matrix in a multiple linear regression.

Examples

### true error term covariance matrix
(Sigma_true <- matrix(c(1,0.5,0.2,0.5,1,0.2,0.2,0.2,2), ncol=3))
### coefficient vector
beta <- matrix(c(-1,1), ncol=1)
### draw data
N <- 100
X <- replicate(N, matrix(rnorm(6), ncol=2), simplify = FALSE)
eps <- replicate(N, rmvnorm(mu = c(0,0,0), Sigma = Sigma_true), simplify = FALSE)
U <- mapply(function(X, eps) X %*% beta + eps, X, eps, SIMPLIFY = FALSE)
### prior parameters for covariance matrix
kappa <- 4
E <- diag(3)
### draw from posterior of coefficient vector
outer_prod <- function(X, U) (U - X %*% beta) %*% t(U - X %*% beta)
S <- Reduce(`+`, mapply(outer_prod, X, U, SIMPLIFY = FALSE))
Sigma_draws <- replicate(100, update_Sigma(kappa, E, N, S))
apply(Sigma_draws, 1:2, mean)
apply(Sigma_draws, 1:2, stats::sd)

RprobitB documentation built on Nov. 10, 2022, 5:12 p.m.