} ) \sim GInvGauss \bigg( 1-k, 1, 2 \sum_{h = 1}^k \frac{|\lambda{jh}|}{\phi_{jh}} \bigg)\end{aligned}$$\
*Step 8* In order to update $\phi_{jh}$ for $j = 1, \cdots, p$ and $h = 1 \cdots, k$, sample $T_{jh}$ from independent Generalized Inverse Gaussian distributions: $$\begin{aligned}
\pi (T_{jh} | \hbox{ --- } ) \sim GInvGauss \bigg( a- 1 , 1 , 2 |\lambda_{jh}| \bigg)\end{aligned}$$\
Then set $\phi_{jh} = \frac{T_{jh}}{\sum_{h = 1}^k T_{jh} }$ *Step 9* Sample $\sigma_j^{-2}$ for $j = 1,\cdots, p$ from conditionally independent gamma distributions $$\begin{aligned}
\pi (\sigma_j^{-2} | \hbox{ --- } ) \sim Gamma \bigg( \frac{1 + n}{2}, \frac{1}{2} + \frac{1}{2} \sum_{i = 1}^n (X_{ij} - \lambda_j^T \eta_i )\bigg)\end{aligned}$$\
]{}

Higher Order Interactions
-------------------------

FIN can be generalized to allow for higher order interactions. In particular, we can obtain estimates for the interaction coefficients up to the $Q^{th}$ order with the following model: $$\begin{aligned}
{\mathbb{E}}(y_i| \eta_i) = \sum_{h= 1}^k \omega_h^{(1)} \eta_{i h} + \sum_{h = 1}^k \omega_h^{(2)} \eta_{i h} ^2 + \cdots + \sum_{h= 1}^k \omega_h^{(Q)} \eta_{i h} ^Q,\end{aligned}$$ which is a polynomial regression in the latent variables. We do not include interactions between the factors, so that the number of parameters to be estimated is $Q  k$. When $Q=2$, this model 