Update class means — update

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

Usage

update_b(beta, Omega, z, m, xi, Dinv)

Arguments

beta: The matrix of the decision-maker specific coefficient vectors of dimension P_r x N. Set to NA if P_r = 0.
Omega: The matrix of class covariance matrices as columns of dimension P_r*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.
xi: The mean vector of length P_r of the normal prior for each b_c. Per default, xi = numeric(P_r).
Dinv: The precision matrix (i.e. the inverse of the covariance matrix) of dimension P_r x P_r of the normal prior for each b_c.

Value

A matrix of updated means for each class in columns.

Details

The following holds independently for each class $c$. Let $b_c$ be the mean of class number $c$. A priori, we assume that $b_c$ is normally distributed with mean vector $\xi$ and covariance matrix $D$. Let $(\beta_n)_{z_n=c}$ be the collection of $\beta_n$ that are currently allocated to class $c$, $m_c$ the class size, and $\bar{b}_c$ their arithmetic mean. Assuming independence across draws, $(\beta_n)_{z_n=c}$ has a normal likelihood of $$\prod_n \phi(\beta_n \mid b_c,\Omega_c),$$ where the product is over the values $n$ for which $z_n=c$ holds. Due to the conjugacy of the prior, the posterior $\Pr(b_c \mid (\beta_n)_{z_n=c})$ follows a normal distribution with mean $$(D^{-1} + m_c\Omega_c^{-1})^{-1}(D^{-1}\xi + m_c\Omega_c^{-1}\bar{b}_c)$$ and covariance matrix $$(D^{-1} + m_c \Omega_c^{-1})^{-1}.$$

Examples

### N = 100 decider, P_r = 2 random coefficients, and C = 2 latent classes
N <- 100
(b_true <- cbind(c(0,0),c(1,1)))
#>      [,1] [,2]
#> [1,]    0    1
#> [2,]    0    1
Omega <- matrix(c(1,0.3,0.3,0.5,1,-0.3,-0.3,0.8), ncol=2)
z <- c(rep(1,N/2),rep(2,N/2))
m <- as.numeric(table(z))
beta <- sapply(z, function(z) rmvnorm(b_true[,z], matrix(Omega[,z],2,2)))
### prior mean vector and precision matrix (inverse of covariance matrix)
xi <- c(0,0)
Dinv <- diag(2)
### updated class means (in columns)
update_b(beta = beta, Omega = Omega, z = z, m = m, xi = xi, Dinv = Dinv)
#>           [,1]      [,2]
#> [1,] 0.1956642 0.7602511
#> [2,] 0.2746751 1.1514358