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

In loelschlaeger/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 (\epsilon_n)_{n=1,\dots,N}.

Details

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

U_n = X_n\beta + \epsilon_n,

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

\Sigma \sim W(\kappa,E)

with \kappa degrees of freedom and scale matrix E. The posterior for \Sigma is the Inverted Wishart distribution with \kappa + N degrees of freedom and scale matrix E^{-1}+S, where S = \sum_{n=1}^{N} \epsilon_n \epsilon_n' is the sum over the outer products of the residuals (\epsilon_n = U_n - X_n\beta)_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)

loelschlaeger/RprobitB documentation built on Oct. 15, 2024, 11:08 a.m.