Update class covariances

This function updates the class covariances (independent from the other classes).

Usage

update_Omega(beta, b, z, m, nu, Theta)

Arguments

beta

The matrix of the decision-maker specific coefficient vectors of dimension P_r x N. Set to NA if P_r = 0.

b

The matrix of class means as columns of dimension P_r x C. Set to NA if P_r = 0.

z

The vector of the allocation variables of length N. Set to NA if P_r = 0.

m

The vector of class sizes of length C.

nu

The degrees of freedom (a natural number greater than P_r) of the Inverse Wishart prior for each Omega_c. Per default, nu = P_r + 2.

Theta

The scale matrix of dimension P_r x P_r of the Inverse Wishart prior for each Omega_c. Per default, Theta = diag(P_r).

Value

A matrix of updated covariance matrices for each class in columns.

Details

The following holds independently for each class \(c\). Let \(\Omega_c\) be the covariance matrix of class number c. A priori, we assume that \(\Omega_c\) is inverse Wishart distributed with \(\nu\) degrees of freedom and scale matrix \(\Theta\). Let \((\beta_n)_{z_n=c}\) be the collection of \(\beta_n\) that are currently allocated to class \(c\), \(m_c\) the size of class \(c\), and \(b_c\) the class mean vector. Due to the conjugacy of the prior, the posterior \(\Pr(\Omega_c \mid (\beta_n)_{z_n=c})\) follows an inverted Wishart distribution with \(\nu + m_c\) degrees of freedom and scale matrix \(\Theta^{-1} + \sum_n (\beta_n - b_c)(\beta_n - b_c)'\), where the product is over the values \(n\) for which \(z_n=c\) holds.

Examples

### N = 100 decider, P_r = 2 random coefficients, and C = 2 latent classes
N <- 100
b <- cbind(c(0,0),c(1,1))
(Omega_true <- matrix(c(1,0.3,0.3,0.5,1,-0.3,-0.3,0.8), ncol=2))
#>      [,1] [,2]
#> [1,]  1.0  1.0
#> [2,]  0.3 -0.3
#> [3,]  0.3 -0.3
#> [4,]  0.5  0.8
z <- c(rep(1,N/2),rep(2,N/2))
m <- as.numeric(table(z))
beta <- sapply(z, function(z) rmvnorm(b[,z], matrix(Omega_true[,z],2,2)))
### degrees of freedom and scale matrix for the Wishart prior
nu <- 1
Theta <- diag(2)
### updated class covariance matrices (in columns)
update_Omega(beta = beta, b = b, z = z, m = m, nu = nu, Theta = Theta)
#>           [,1]       [,2]
#> [1,] 0.6972189  1.0776332
#> [2,] 0.2217133 -0.2271837
#> [3,] 0.2217133 -0.2271837
#> [4,] 0.3307745  1.3250835