# Distribution of factor scores

Given that $X_i\sim N_p(\mu,LL^T+\Psi)$ where $\mu_{p\times p}$ is the mean vector $L_{p\times p}$ is the loading matrix and $\Psi_{p\times p}$ is the matrix of specific variances. Find the joint distribution of $(F_i^M;F_i^R)$ where $$F_i^M=(L^T\Psi^{-1}L)^{-1}L^T\Psi^{-1}(X_i-\mu)$$ $$F_i^R=L^T(LL^T+\Psi)^{-1}(X_i-\mu)$$

Following (8.14) of Matrix Cookbook

$$F_i^M\sim N_p(0_p,(L^T\Psi^{-1}L)^{-1}L^T\Psi^{-1}(LL^T+\Psi)((L^T\Psi^{-1}L)^{-1}L^T\Psi^{-1})^T)$$ $$F_i^R\sim N_p(0,L^T(LL^T+\Psi)^{-1}(LL^T+\Psi)(L^T(LL^T+\Psi)^{-1})^T)$$

Is this right? How I can find the joint distribution?

Yes, the distribution of $F_i^M$ and $F_i^R$ are correct because if $x \sim \mathcal{N}(\mu_n, \Sigma_{n \times n})$, then the distribution of $A x$, where $A \in \mathbb{R}^{n \times n}$ is of full rank, is given as $Ax \sim \mathcal{N}(A\mu, A\Sigma A^T)$.
Notice that both $F_i^M$ and $F_i^R$ are dependent on $X_i$, i.e., we can write the distribution of \begin{align} F_i^M(X_i) &\sim \mathcal{N}(0_p, Q \Gamma Q^T) \\ F_i^R(X_i) &\sim \mathcal{N}(0_p, R \Gamma R^T), \\ \end{align} where $\Gamma = L L^T + \Psi$, $Q = (L^T \Psi^{-1} L)^{-1}L^T \Psi^{-1}$ and $R = L^T(L L^T + \Psi)^{-1}$. As, these two are not independent random variables, we need to find $\mathbb{E}(F_i^M (F_i^R)^T)$ and $\mathbb{E}(F_i^R (F_i^M)^T)$.
In the same way, $\mathbb{E}(F_i^R (F_i^M)^T) = R \Gamma Q^T$. Now, combining all the terms, we get the joint distribution as follows: