Via the normalization constant we have

$$
\begin{align*}
& \int \mathcal{N}_L (\hat{\mathbf{x}} | \hat{\mathbf{m}}_k, (1 + \hat{\lambda}_k^{-1}) \Gamma_k^{-1}, \mathcal{X}) W (\Gamma_k | \hat{\Psi}_k, \hat{\nu}_k, \mathcal{X}) d\Gamma_k \\
&= \frac{C(\hat{\Psi}_k, \hat{\nu}_k, \hat{\lambda}_k)}{B(\hat{\Psi}_k^*, \hat{\nu}_k^*)} \underbrace{\int W (\Gamma_k | \Psi_k^*, \hat{\nu}_k^*) d\Gamma_k}_{=1} = \frac{C(\hat{\Psi}_k, \hat{\nu}_k, \hat{\lambda}_k)}{B(\hat{\Psi}_k^*, \hat{\nu}_k^*)} = \frac{B(\hat{\Psi}_k, \hat{\nu}_k)}{(2\pi (1 + \hat{\lambda}_k^{-1}))^{L/2}} \frac{1}{B(\hat{\Psi}_k^*, \hat{\nu}_k^*)} \\
&= \frac{1}{(2\pi (1 + \hat{\lambda}_k^{-1}))^{L/2}} \frac{|\hat{\Psi}_k|^{-\hat{\nu}_k/2} (2^{\hat{\nu}_k L/2} \pi^{L(L-1)/4} \prod_{l=1}^L \Gamma(\frac{\hat{\nu}_k+1-l}{2}))^{-1}}{|\hat{\Psi}_k^*|^{-(\hat{\nu}_k+1)/2} (2^{(\hat{\nu}_k+1)L/2} \pi^{L(L-1)/4} \prod_{l=1}^L \Gamma(\frac{(\hat{\nu}_k+1)+1-l}{2}))^{-1}} \\
&= \frac{1}{(\pi (1 + \hat{\lambda}_k^{-1}))^{L/2}} \frac{|\hat{\Psi}_k|^{-\hat{\nu}_k/2} \Gamma(\frac{\hat{\nu}_k+1}{2}) \Gamma(\frac{\hat{\nu}_k}{2}) \dots \Gamma(\frac{\hat{\nu}_k+2-L}{2})}{|\hat{\Psi}_k^*|^{-(\hat{\nu}_k+1)/2} \Gamma(\frac{\hat{\nu}_k}{2}) \Gamma(\frac{\hat{\nu}_k-1}{2}) \dots \Gamma(\frac{\hat{\nu}_k+2-L}{2}) \Gamma(\frac{\hat{\nu}_k+1-L}{2})} \\
&= \frac{\Gamma(\frac{\hat{\nu}_k+1}{2})}{\Gamma(\frac{\hat{\nu}_k+1-L}{2}) \pi^{L/2}} \frac{|\hat{\Psi}_k|^{-\hat{\nu}_k/2}}{(1 + \hat{\lambda}_k^{-1})^{L/2}} |\hat{\Psi}_k|^{(\hat{\nu}_k+1)/2} [1 + (1 + \hat{\lambda}_k^{-1})^{-1} (\hat{\mathbf{x}} - \hat{\mathbf{m}}_k)^{\top} \hat{\Psi}_k (\hat{\mathbf{x}} - \hat{\mathbf{m}}_k)]^{-(\hat{\nu}_k+1)/2} \\
&= \frac{\Gamma(\frac{\hat{\nu}_k+1}{2})}{\Gamma(\frac{\hat{\nu}_k+1-L}{2}) \pi^{L/2}} \frac{|\hat{\Psi}_k|^{1/2}}{(1 + \hat{\lambda}_k^{-1})^{L/2}} [1 + (1 + \hat{\lambda}_k^{-1})^{-1} (\hat{\mathbf{x}} - \hat{\mathbf{m}}_k)^{\top} \hat{\Psi}_k (\hat{\mathbf{x}} - \hat{\mathbf{m}}_k)]^{-(\hat{\nu}_k+1)/2} \\
&= \text{St}(\hat{\mathbf{x}} | \hat{\mathbf{m}}_k, \mathbf{L}_k, \hat{\nu}_k + 1 - L).
\end{align*}
$$

Here,

$$
\mathbf{L}_k = \frac{(\hat{\nu}_k + 1 - L) \hat{\lambda}_k}{1 + \hat{\lambda}_k} \hat{\Psi}_k,
$$

and $\Delta^2$ is the squared Mahalanobis distance defined by

$$
\Delta^2 = (\hat{\mathbf{x}} - \hat{\mathbf{m}}_k)^{\top} \mathbf{L}_k (\hat{\mathbf{x}} - \hat{\mathbf{m}}_k).
$$
Then, the last equality holds since we have

$$
\begin{align*}
\mathrm{St}(\hat{\mathbf{x}} | \hat{\mathbf{m}}_k, \mathbf{L}_k, \hat{\nu}_k + 1 - L) &= \frac{\Gamma\left(\frac{\hat{\nu}_k + 1 - L}{2} + \frac{L}{2}\right)}{\Gamma\left(\frac{\hat{\nu}_k + 1 - L}{2}\right) \pi^{L/2}} \frac{|\mathbf{L}_k|^{1/2}}{(\hat{\nu}_k + 1 - L)^{L/2}} \left[ 1 + \frac{\Delta^2}{\hat{\nu}_k + 1 - L} \right]^{-\frac{(\hat{\nu}_k + 1 - L)/2 - L/2}{2}} \\
&= \frac{\Gamma\left(\frac{\hat{\nu}_k + 1}{2}\right)}{\Gamma\left(\frac{\hat{\nu}_k + 1 - L}{2}\right) \pi^{L/2}} \frac{(\hat{\nu}_k + 1 - L)^{L/2} |\hat{\Psi}_k|^{1/2}}{(\hat{\nu}_k + 1 - L)^{L/2} (1 + \hat{\lambda}_k^{-1})^{L/2}} \left[ 1 + \frac{\Delta^2}{\hat{\nu}_k + 1 - L} \right]^{-\frac{(\hat{\nu}_k + 1 - L)/2 - L/2}{2}} \\
&= \frac{\Gamma\left(\frac{\hat{\nu}_k + 1}{2}\right)}{\Gamma\left(\frac{\hat{\nu}_k + 1 - L}{2}\right) \pi^{L/2}} \frac{|\hat{\Psi}_k|^{1/2}}{(1 + \hat{\lambda}_k^{-1})^{L/2}} \left[ 1 + \frac{(\hat{\nu}_k + 1 - L) \hat{\lambda}_k (\hat{\mathbf{x}} - \hat{\mathbf{m}}_k)^{\top} \hat{\Psi}_k (\hat{\mathbf{x}} - \hat{\mathbf{m}}_k)}{(1 + \hat{\lambda}_k) \hat{\nu}_k + 1 - L} \right]^{-\frac{(\hat{\nu}_k + 1)}{2}} \\
&= \frac{\Gamma\left(\frac{\hat{\nu}_k + 1}{2}\right)}{\Gamma\left(\frac{\hat{\nu}_k + 1 - L}{2}\right) \pi^{L/2}} \frac{|\hat{\Psi}_k|^{1/2}}{(1 + \hat{\lambda}_k^{-1})^{L/2}} \left[ 1 + \frac{(\hat{\mathbf{x}} - \hat{\mathbf{m}}_k)^{\top} \hat{\Psi}_k (\hat{\mathbf{x}} - \hat{\mathbf{m}}_k)}{(1 + \hat{\lambda}_k^{-1})} \right]^{-\frac{(\hat{\nu}_k + 1)}{2}}.
\end{align*}
$$

F.5. *Proof of Theorem F.3*

From the product rule of probability, we see that this conditional distribution can be
evaluated from the joint and marginal distributions. Furthermore, by integrating out $\hat{z}$
and $\Theta$, the predictive conditional density is then given by

$$
p(\hat{\mathbf{y}} | \hat{\mathbf{x}}, \mathcal{X}, \mathcal{Y}) = \frac{p(\hat{\mathbf{y}}, \hat{\mathbf{x}} | \mathcal{X}, \mathcal{Y})}{p(\hat{\mathbf{x}} | \mathcal{X}, \mathcal{Y})} = \frac{\sum_{\hat{z}} \int p(\hat{\mathbf{y}}, \hat{\mathbf{x}}, \hat{z}, \Theta | \mathcal{X}, \mathcal{Y}) d\hat{z} d\Theta}{\sum_{\hat{z}} \int p(\hat{\mathbf{x}}, \hat{z}, \Theta | \mathcal{X}, \mathcal{Y}) d\hat{z} d\Theta} \equiv \frac{T_1}{T_2}.
$$

Next, with a similar step as in the proof of Theorem F.2, we also obtain

$$
T_2 = \sum_{k=1}^{K} \mathbb{E}_{q_{\tau}} [\pi_k(\tau)] \mathrm{St}(\hat{\mathbf{x}} | \hat{\mathbf{m}}_k, \mathbf{L}_k, \hat{\nu}_k + 1 - L).
$$

Therefore

$$
\begin{align*}
p(\hat{\mathbf{y}} | \hat{\mathbf{x}}, \mathcal{X}, \mathcal{Y})
  &\approx \frac{\sum_{k=1}^{K} \mathbb{E}_{q_{\tau}} [\pi_k(\tau)] \mathrm{St}(\hat{\mathbf{x}} | \hat{\mathbf{m}}_k, \mathbf{L}_k, \hat{\nu}_k + 1 - L) \mathcal{N}_D(\hat{\mathbf{y}} | \hat{\mathbf{A}}_k \hat{\mathbf{x}} + \hat{\mathbf{b}}_k, \hat{\mathbf{\Sigma}}_k)}{\sum_{k=1}^{K} \mathbb{E}_{q_{\tau}} [\pi_k(\tau)] \mathrm{St}(\hat{\mathbf{x}} | \hat{\mathbf{m}}_k, \mathbf{L}_k, \hat{\nu}_k + 1 - L)} \\
  &= \sum_{k=1}^{K} \frac{\mathbb{E}_{q_{\tau}} [\pi_k(\tau)] \mathrm{St}(\hat{\mathbf{x}} | \hat{\mathbf{m}}_k, \mathbf{L}_k, \hat{\nu}_k + 1 - L)}{\sum_{k=1}^{K} \mathbb{E}_{q_{\tau}} [\pi_k(\tau)] \mathrm{St}(\hat{\mathbf{x}} | \hat{\mathbf{m}}_k, \mathbf{L}_k, \hat{\nu}_k + 1 - L)} \mathcal{N}_D(\hat{\mathbf{y}} | \hat{\mathbf{A}}_k \hat{\mathbf{x}} + \hat{\mathbf{b}}_k, \hat{\mathbf{\Sigma}}_k) \\
  &\equiv \sum_{k=1}^{K} g_k(\hat{\mathbf{x}} | \hat{\Theta}, \hat{\phi}, \mathcal{X}, \mathcal{Y}) \mathcal{N}_D(\hat{\mathbf{y}} | \hat{\mathbf{A}}_k \hat{\mathbf{x}} + \hat{\mathbf{b}}_k, \hat{\mathbf{\Sigma}}_k),
\end{align*}
$$

which is a mixture of Gaussian experts since we have

$$
g_k(\hat{\mathbf{x}} | \hat{\Theta}, \hat{\phi}, \mathcal{X}, \mathcal{Y}) = \frac{\mathbb{E}_{q_{\tau}}[\pi_k(\tau)] \mathrm{St}(\hat{\mathbf{x}} | \hat{\mathbf{m}}_k, \mathbf{L}_k, \hat{\nu}_k + 1 - L)}{\sum_{k=1}^{K} \mathbb{E}_{q_{\tau}}[\pi_k(\tau)] \mathrm{St}(\hat{\mathbf{x}} | \hat{\mathbf{m}}_k, \mathbf{L}_k, \hat{\nu}_k + 1 - L)}, \quad k \in [K],
$$

belongs to a $K - 1$ dimensional probability simplex.
Then, the last equality holds since we have

$$
\begin{align*}
\mathrm{St}(\hat{\mathbf{x}} | \hat{\mathbf{m}}_k, \mathbf{L}_k, \hat{\nu}_k + 1 - L) &= \frac{\Gamma\left(\frac{\hat{\nu}_k+1-L}{2} + \frac{L}{2}\right)}{\Gamma\left(\frac{\hat{\nu}_k+1-L}{2}\right) \pi^{L/2}} \frac{|\mathbf{L}_k|^{1/2}}{(\hat{\nu}_k + 1 - L)^{L/2}} \left[1 + \frac{\Delta^2}{\hat{\nu}_k + 1 - L}\right]^{-\frac{(\hat{\nu}_k+1-L)}{2}-L/2} \\
&= \frac{\Gamma\left(\frac{\hat{\nu}_k+1}{2}\right)}{\Gamma\left(\frac{\hat{\nu}_k+1-L}{2}\right) \pi^{L/2}} \frac{(\hat{\nu}_k + 1 - L)^{L/2} |\hat{\Psi}_k|^{1/2}}{(\hat{\nu}_k + 1 - L)^{L/2} (1 + \hat{\lambda}_k^{-1})^{L/2}} \left[1 + \frac{\Delta^2}{\hat{\nu}_k + 1 - L}\right]^{-\frac{(\hat{\nu}_k+1-L)}{2}-L/2} \\
&= \frac{\Gamma\left(\frac{\hat{\nu}_k+1}{2}\right)}{\Gamma\left(\frac{\hat{\nu}_k+1-L}{2}\right) \pi^{L/2}} \frac{|\hat{\Psi}_k|^{1/2}}{(1 + \hat{\lambda}_k^{-1})^{L/2}} \left[1 + \frac{(\hat{\nu}_k + 1 - L) \hat{\lambda}_k (\hat{\mathbf{x}} - \hat{\mathbf{m}}_k)^{\top} \hat{\Psi}_k (\hat{\mathbf{x}} - \hat{\mathbf{m}}_k)}{(1 + \hat{\lambda}_k) \hat{\nu}_k + 1 - L}\right]^{-\frac{(\hat{\nu}_k+1)}{2}} \\
&= \frac{\Gamma\left(\frac{\hat{\nu}_k+1}{2}\right)}{\Gamma\left(\frac{\hat{\nu}_k+1-L}{2}\right) \pi^{L/2}} \frac{|\hat{\Psi}_k|^{1/2}}{(1 + \hat{\lambda}_k^{-1})^{L/2}} \left[1 + \frac{(\hat{\mathbf{x}} - \hat{\mathbf{m}}_k)^{\top} \hat{\Psi}_k (\hat{\mathbf{x}} - \hat{\mathbf{m}}_k)}{(1 + \hat{\lambda}_k^{-1})}\right]^{-\frac{(\hat{\nu}_k+1)}{2}}.
\end{align*}
$$

F.5. *Proof of Theorem F.3*

From the product rule of probability, we see that this conditional distribution can be
evaluated from the joint and marginal distributions. Furthermore, by integrating out $\hat{z}$
and $\Theta$, the predictive conditional density is then given by

$$
p(\hat{\mathbf{y}} | \hat{\mathbf{x}}, \mathcal{X}, \mathcal{Y}) = \frac{p(\hat{\mathbf{y}}, \hat{\mathbf{x}} | \mathcal{X}, \mathcal{Y})}{p(\hat{\mathbf{x}} | \mathcal{X}, \mathcal{Y})} = \frac{\sum_{\hat{z}} \int p(\hat{\mathbf{y}}, \hat{\mathbf{x}}, \hat{z}, \Theta | \mathcal{X}, \mathcal{Y}) d\hat{z} d\Theta}{\sum_{\hat{z}} \int p(\hat{\mathbf{x}}, \hat{z}, \Theta | \mathcal{X}, \mathcal{Y}) d\hat{z} d\Theta} \equiv \frac{T_1}{T_2}.
$$

Next, with a similar step as in the proof of Theorem F.2, we also obtain

$$
T_2 = \sum_{k=1}^{K} \mathbb{E}_{q_{\tau}} [\pi_k(\tau)] \mathrm{St}(\hat{\mathbf{x}} | \hat{\mathbf{m}}_k, \mathbf{L}_k, \hat{\nu}_k + 1 - L).
$$

Therefore

$$
\begin{align*}
p(\hat{\mathbf{y}} | \hat{\mathbf{x}}, \mathcal{X}, \mathcal{Y})
  &\approx \frac{\sum_{k=1}^{K} \mathbb{E}_{q_{\tau}} [\pi_k(\tau)] \mathrm{St}(\hat{\mathbf{x}} | \hat{\mathbf{m}}_k, \mathbf{L}_k, \hat{\nu}_k + 1 - L) \mathcal{N}_D(\hat{\mathbf{y}} | \hat{\mathbf{A}}_k \hat{\mathbf{x}} + \hat{\mathbf{b}}_k, \hat{\mathbf{\Sigma}}_k)}{\sum_{k=1}^{K} \mathbb{E}_{q_{\tau}} [\pi_k(\tau)] \mathrm{St}(\hat{\mathbf{x}} | \hat{\mathbf{m}}_k, \mathbf{L}_k, \hat{\nu}_k + 1 - L)} \\
  &= \sum_{k=1}^{K} \frac{\mathbb{E}_{q_{\tau}} [\pi_k(\tau)] \mathrm{St}(\hat{\mathbf{x}} | \hat{\mathbf{m}}_k, \mathbf{L}_k, \hat{\nu}_k + 1 - L)}{\sum_{k=1}^{K} \mathbb{E}_{q_{\tau}} [\pi_k(\tau)] \mathrm{St}(\hat{\mathbf{x}} | \hat{\mathbf{m}}_k, \mathbf{L}_k, \hat{\nu}_k + 1 - L)} \mathcal{N}_D(\hat{\mathbf{y}} | \hat{\mathbf{A}}_k \hat{\mathbf{x}} + \hat{\mathbf{b}}_k, \hat{\mathbf{\Sigma}}_k) \\
  &\equiv \sum_{k=1}^{K} g_k(\hat{\mathbf{x}} | \hat{\Theta}, \hat{\phi}, \mathcal{X}, \mathcal{Y}) \mathcal{N}_D(\hat{\mathbf{y}} | \hat{\mathbf{A}}_k \hat{\mathbf{x}} + \hat{\mathbf{b}}_k, \hat{\mathbf{\Sigma}}_k),
\end{align*}
$$

which is a mixture of Gaussian experts since we have

$$
g_k(\hat{\mathbf{x}} | \hat{\Theta}, \hat{\phi}, \mathcal{X}, \mathcal{Y}) = \frac{\mathbb{E}_{q_{\tau}}[\pi_k(\tau)] \mathrm{St}(\hat{\mathbf{x}} | \hat{\mathbf{m}}_k, \mathbf{L}_k, \hat{\nu}_k + 1 - L)}{\sum_{k=1}^{K} \mathbb{E}_{q_{\tau}}[\pi_k(\tau)] \mathrm{St}(\hat{\mathbf{x}} | \hat{\mathbf{m}}_k, \mathbf{L}_k, \hat{\nu}_k + 1 - L)}, \quad k \in [K],
$$

belongs to a $K-1$ dimensional probability simplex.
### F.6. Proof of Theorem F.4

To deal with high-dimensional regression data, namely high-to-low regression, given the inverse conditional density $p(\hat{y} | \hat{x}, \mathcal{X}, \mathcal{Y})$, we want to compute the following forward conditional density

$$p(\hat{x} | \hat{y}, \mathcal{X}, \mathcal{Y}) = \frac{p(\hat{x}, \hat{y} | \mathcal{X}, \mathcal{Y})}{p(\hat{y} | \mathcal{X}, \mathcal{Y})} = \frac{p(\hat{x}, \hat{y} | \mathcal{X}, \mathcal{Y})}{\int_{\hat{x}} p(\hat{x}, \hat{y} | \mathcal{X}, \mathcal{Y}) d\hat{x}} = \frac{T_1}{\int_{\hat{x}} T_1(\hat{x}) d\hat{x}} \equiv \frac{T_1}{T_3}.$$

Then, we have to compute or numerically approximate $D_3$. Using Theorem F.2 and definition of Student's t-distribution, we obtain

$$T_3 = \sum_{k=1}^{K} \mathbb{E}_{q_{\tau}} [\pi_k(\tau)] D_k.$$

Then, by definition of Student's t-distribution, it holds that

$$
\begin{align*}
D_k &= \int \mathrm{St}(\hat{x} | \hat{\mathbf{m}}_k, \mathbf{L}_k, \hat{\nu}_k + 1 - L) \mathcal{N}_D (\hat{y} | \hat{\mathbf{A}}_k \hat{x} + \hat{\mathbf{b}}_k, \hat{\mathbf{\Sigma}}_k) d\hat{x} \\
&= \int \int_0^\infty \mathcal{N}_L (\hat{x} | \hat{\mathbf{m}}_k, (\eta \mathbf{L}_k)^{-1}) \mathrm{Gam} \left( \eta \left| \frac{\hat{\nu}_k + 1 - L}{2} \right., \frac{\hat{\nu}_k + 1 - L}{2} \right) \mathcal{N}_D (\hat{y} | \hat{\mathbf{A}}_k \hat{x} + \hat{\mathbf{b}}_k, \hat{\mathbf{\Sigma}}_k) d\eta d\hat{x} \\
&= \int_0^\infty \int \mathcal{N}_D (\hat{y} | \hat{\mathbf{A}}_k \hat{x} + \hat{\mathbf{b}}_k, \hat{\mathbf{\Sigma}}_k) \mathcal{N}_L (\hat{x} | \hat{\mathbf{m}}_k, (\eta \mathbf{L}_k)^{-1}) d\hat{x} \mathrm{Gam} \left( \eta \left| \frac{\hat{\nu}_k + 1 - L}{2} \right., \frac{\hat{\nu}_k + 1 - L}{2} \right) d\eta \\
&= \int_0^\infty \mathcal{N}_D (\hat{y} | \hat{\mathbf{A}}_k \hat{\mathbf{m}}_k + \hat{\mathbf{b}}_k, \hat{\mathbf{\Sigma}}_k + \eta^{-1} \hat{\mathbf{A}}_k \mathbf{L}_k^{-1} \hat{\mathbf{A}}_k^\top) \mathrm{Gam} \left( \eta \left| \frac{\hat{\nu}_k + 1 - L}{2} \right., \frac{\hat{\nu}_k + 1 - L}{2} \right) d\eta.
\end{align*}
$$

Furthermore, we used the fact that

$$\int \mathcal{N}_D (\hat{y} | \hat{\mathbf{A}}_k \hat{x} + \hat{\mathbf{b}}_k, \hat{\mathbf{\Sigma}}_k) \mathcal{N}_L (\hat{x} | \hat{\mathbf{m}}_k, (\eta \mathbf{L}_k)^{-1}) d\hat{x} = \mathcal{N}_D (\hat{y} | \hat{\mathbf{A}}_k \hat{\mathbf{m}}_k + \hat{\mathbf{b}}_k, \hat{\mathbf{\Sigma}}_k + \hat{\mathbf{A}}_k (\eta \mathbf{L}_k)^{-1} \hat{\mathbf{A}}_k^\top).$$

Indeed, we made use of the following results for marginal and conditional Gaussians, see, e.g., Bishop (2006, Eq. (2.115), page 93). Given a marginal Gaussian distribution for x and a conditional Gaussian distribution for y given x in the form

$$\begin{align*}
p(\mathbf{x}) &= \mathcal{N}(\mathbf{x} | \boldsymbol{\mu}, \boldsymbol{\Gamma}^{-1}), \\
p(\mathbf{y} | \mathbf{x}) &= \mathcal{N}(\mathbf{y} | \mathbf{A}\mathbf{x} + \mathbf{b}, \mathbf{L}^{-1}),
\end{align*}$$

then the marginal distribution of y and the conditional distribution of x given y are given by

$$\begin{align*}
p(y) &= \int p(y | x)p(x)dx = \mathcal{N}(y | \mathbf{A}\boldsymbol{\mu} + \mathbf{b}, \mathbf{L}^{-1} + \mathbf{A}\boldsymbol{\Gamma}^{-1}\mathbf{A}^{\top}), \\
p(x | y) &= \mathcal{N}(x | \boldsymbol{\Sigma}\{\mathbf{A}^{\top}\mathbf{L}(y - \mathbf{b}) + \boldsymbol{\Gamma}\boldsymbol{\mu}\}, \boldsymbol{\Sigma}), \quad \boldsymbol{\Sigma} = (\boldsymbol{\Gamma} + \mathbf{A}^{\top}\mathbf{L}\mathbf{A})^{-1}.
\end{align*}$$

In our situation, the desired result is obtained via using $y \equiv \hat{y}, \mathbf{A} \equiv \hat{\mathbf{A}}_k, \mathbf{b} \equiv \hat{\mathbf{b}}_k, \mathbf{L}^{-1} = \hat{\mathbf{\Sigma}}_k, x \equiv \hat{x}, \boldsymbol{\mu} \equiv \hat{\mathbf{m}}_k, \boldsymbol{\Gamma}^{-1} \equiv (\eta \mathbf{L}_k)^{-1}$.
F.6. Proof of Theorem F.4

To deal with high-dimensional regression data, namely high-to-low regression, given the inverse conditional density $p(\hat{y} | \hat{x}, \mathcal{X}, \mathcal{Y})$, we want to compute the following forward conditional density

$$p(\hat{x} | \hat{y}, \mathcal{X}, \mathcal{Y}) = \frac{p(\hat{x}, \hat{y} | \mathcal{X}, \mathcal{Y})}{p(\hat{y} | \mathcal{X}, \mathcal{Y})} = \frac{p(\hat{x}, \hat{y} | \mathcal{X}, \mathcal{Y})}{\int_{\hat{x}} p(\hat{x}, \hat{y} | \mathcal{X}, \mathcal{Y}) d\hat{x}} = \frac{T_1}{\int_{\hat{x}} T_1(\hat{x}) d\hat{x}} \equiv \frac{T_1}{T_3}.$$

Then, we have to compute or numerically approximate $D_3$. Using Theorem F.2 and definition of Student's t-distribution, we obtain

$$T_3 = \sum_{k=1}^{K} \mathbb{E}_{q_{\tau}} [\pi_k(\tau)] D_k.$$

Then, by definition of Student's t-distribution, it holds that

$$\begin{align*}
D_k &= \int \mathrm{St}(\hat{x} | \hat{\mathbf{m}}_k, \mathbf{L}_k, \hat{\nu}_k + 1 - L) \mathcal{N}_D (\hat{y} | \hat{\mathbf{A}}_k \hat{x} + \hat{\mathbf{b}}_k, \hat{\mathbf{\Sigma}}_k) d\hat{x} \\
&= \int \int_0^\infty \mathcal{N}_L (\hat{x} | \hat{\mathbf{m}}_k, (\eta \mathbf{L}_k)^{-1}) \mathrm{Gam} \left( \eta \left| \frac{\hat{\nu}_k + 1 - L}{2} \right., \frac{\hat{\nu}_k + 1 - L}{2} \right) \mathcal{N}_D (\hat{y} | \hat{\mathbf{A}}_k \hat{x} + \hat{\mathbf{b}}_k, \hat{\mathbf{\Sigma}}_k) d\eta d\hat{x} \\
&= \int_0^\infty \int \mathcal{N}_D (\hat{y} | \hat{\mathbf{A}}_k \hat{x} + \hat{\mathbf{b}}_k, \hat{\mathbf{\Sigma}}_k) \mathcal{N}_L (\hat{x} | \hat{\mathbf{m}}_k, (\eta \mathbf{L}_k)^{-1}) d\hat{x} \mathrm{Gam} \left( \eta \left| \frac{\hat{\nu}_k + 1 - L}{2} \right., \frac{\hat{\nu}_k + 1 - L}{2} \right) d\eta \\
&= \int_0^\infty \mathcal{N}_D (\hat{y} | \hat{\mathbf{A}}_k \hat{\mathbf{m}}_k + \hat{\mathbf{b}}_k, \hat{\mathbf{\Sigma}}_k + \eta^{-1} \hat{\mathbf{A}}_k \mathbf{L}_k^{-1} \hat{\mathbf{A}}_k^\top) \mathrm{Gam} \left( \eta \left| \frac{\hat{\nu}_k + 1 - L}{2} \right., \frac{\hat{\nu}_k + 1 - L}{2} \right) d\eta.
\end{align*}$$

Furthermore, we used the fact that

$$\int \mathcal{N}_D (\hat{y} | \hat{\mathbf{A}}_k \hat{x} + \hat{\mathbf{b}}_k, \hat{\mathbf{\Sigma}}_k) \mathcal{N}_L (\hat{x} | \hat{\mathbf{m}}_k, (\eta \mathbf{L}_k)^{-1}) d\hat{x} = \mathcal{N}_D (\hat{y} | \hat{\mathbf{A}}_k \hat{\mathbf{m}}_k + \hat{\mathbf{b}}_k, \hat{\mathbf{\Sigma}}_k + \hat{\mathbf{A}}_k (\eta \mathbf{L}_k)^{-1} \hat{\mathbf{A}}_k^\top).$$

Indeed, we made use of the following results for marginal and conditional Gaussians, see, e.g., Bishop (2006, Eq. (2.115), page 93). Given a marginal Gaussian distribution for x and a conditional Gaussian distribution for y given x in the form

$$\begin{align*}
p(\mathbf{x}) &= \mathcal{N}(\mathbf{x} | \boldsymbol{\mu}, \boldsymbol{\Gamma}^{-1}), \\
p(\mathbf{y} | \mathbf{x}) &= \mathcal{N}(\mathbf{y} | \mathbf{A}\mathbf{x} + \mathbf{b}, \mathbf{L}^{-1}),
\end{align*}$$

then the marginal distribution of y and the conditional distribution of x given y are
given by

$$\begin{align*}
p(\mathbf{y}) &= \int p(\mathbf{y} | \mathbf{x}) p(\mathbf{x}) d\mathbf{x} = \mathcal{N} \left( \mathbf{y} | \mathbf{A}\boldsymbol{\mu} + \mathbf{b}, \mathbf{L}^{-1} + \mathbf{A}\boldsymbol{\Gamma}^{-1}\mathbf{A}^{\top} \right), \\
p(\mathbf{x} | \mathbf{y}) &= \mathcal{N} \left( \mathbf{x} | \boldsymbol{\Sigma} \left\{ \mathbf{A}^{\top}\mathbf{L}(\mathbf{y} - \mathbf{b}) + \boldsymbol{\Gamma}\boldsymbol{\mu} \right\}, \boldsymbol{\Sigma} \right), \boldsymbol{\Sigma} = (\boldsymbol{\Gamma} + \mathbf{A}^{\top}\mathbf{L}\mathbf{A})^{-1}.
\end{align*}$$

In our situation, the desired result is obtained via using $\mathbf{y} \equiv \hat{\mathbf{y}}$, $\mathbf{A} \equiv \hat{\mathbf{A}}_k$, $\mathbf{b} \equiv \hat{\mathbf{b}}_k$, $\mathbf{L}^{-1} = \hat{\mathbf{\Sigma}}_k$, $\mathbf{x} \equiv \hat{\mathbf{x}}$, $\boldsymbol{\mu} \equiv \hat{\mathbf{m}}_k$, $\boldsymbol{\Gamma}^{-1} \equiv (\eta \mathbf{L}_k)^{-1}$.
Therefore, we obtain

$p(\hat{x} | \hat{y}, \mathcal{X}, \mathcal{Y})$

$$
\begin{align*}
&\approx \sum_{k=1}^{K} \frac{\mathbb{E}_{q_{\tau}}[\pi_k(\tau)] \mathrm{St}(\hat{x} | \hat{\mathbf{m}}_k, \mathbf{L}_k, \hat{\nu}_k + 1 - L) \mathcal{N}_D(\hat{y} | \hat{\mathbf{A}}_k \hat{\mathbf{x}} + \hat{\mathbf{b}}_k, \hat{\mathbf{\Sigma}}_k)}{\sum_{k=1}^{K} \mathbb{E}_{q_{\tau}}[\pi_k(\tau)] \int \mathrm{St}(\hat{x} | \hat{\mathbf{m}}_k, \mathbf{L}_k, \hat{\nu}_k + 1 - L) \mathcal{N}_D(\hat{y} | \hat{\mathbf{A}}_k \hat{\mathbf{x}} + \hat{\mathbf{b}}_k, \hat{\mathbf{\Sigma}}_k) d\hat{\mathbf{x}}} \\
&= \sum_{k=1}^{K} \frac{\mathbb{E}_{q_{\tau}}[\pi_k(\tau)] \int_{0}^{\infty} \mathcal{N}_D(\hat{y} | \hat{\mathbf{A}}_k \hat{\mathbf{x}} + \hat{\mathbf{b}}_k, \hat{\mathbf{\Sigma}}_k) \mathcal{N}_L(\hat{x} | \hat{\mathbf{m}}_k, (\eta \mathbf{L}_k)^{-1}) \mathrm{Gam}(\eta | \frac{\hat{\nu}_k + 1 - L}{2}, \frac{\hat{\nu}_k + 1 - L}{2}) d\eta}{\sum_{k=1}^{K} \mathbb{E}_{q_{\tau}}[\pi_k(\tau)] \int_{0}^{\infty} \mathcal{N}_D(\hat{y} | \hat{\mathbf{A}}_k \hat{\mathbf{m}}_k + \hat{\mathbf{b}}_k, \hat{\mathbf{\Sigma}}_k + \eta^{-1} \hat{\mathbf{A}}_k \mathbf{L}_k^{-1} \hat{\mathbf{A}}_k^{\top}) \mathrm{Gam}(\eta | \frac{\hat{\nu}_k + 1 - L}{2}, \frac{\hat{\nu}_k + 1 - L}{2}) d\eta} \\
&= \sum_{k=1}^{K} \frac{\mathbb{E}_{q_{\tau}}[\pi_k(\tau)] \int_{0}^{\infty} \mathcal{N}_D(\hat{y} | \hat{\mathbf{c}}_k^*, \hat{\mathbf{\Gamma}}_k^*(\eta)) \mathcal{N}_L(\hat{x} | \hat{\mathbf{A}}_k^*(\eta)\hat{y} + \hat{\mathbf{b}}_k^*(\eta), \hat{\mathbf{\Sigma}}_k^*(\eta)) \mathrm{Gam}(\eta | \frac{\hat{\nu}_k + 1 - L}{2}, \frac{\hat{\nu}_k + 1 - L}{2}) d\eta}{\sum_{k=1}^{K} \mathbb{E}_{q_{\tau}}[\pi_k(\tau)] \int_{0}^{\infty} \mathcal{N}_D(\hat{y} | \hat{\mathbf{c}}_k^*, \hat{\mathbf{\Gamma}}_k^*(\eta)) \mathrm{Gam}(\eta | \frac{\hat{\nu}_k + 1 - L}{2}, \frac{\hat{\nu}_k + 1 - L}{2}) d\eta} \\
&\approx \sum_{i=1}^{I} \sum_{k=1}^{K} g_{ki}(\hat{y} | \hat{\Theta}^*, \hat{\phi}^*, \mathcal{X}, \mathcal{Y}) \mathcal{N}_L(\hat{x} | \hat{\mathbf{A}}_k^*(\eta_i)\hat{y} + \hat{\mathbf{b}}_k^*(\eta_i), \hat{\mathbf{\Sigma}}_k^*(\eta_i)),
\end{align*}
$$

where, for all $k \in [K], i \in [I]$,

$$
g_{ki}(\hat{y} | \hat{\Theta}^*, \hat{\phi}^*, \mathcal{X}, \mathcal{Y}) = \frac{\mathbb{E}_{q_{\tau}}[\pi_k(\tau)] \mathcal{N}_D(\hat{y} | \hat{\mathbf{c}}_k^*, \hat{\mathbf{\Gamma}}_k^*(\eta_i)) \mathrm{Gam}(\eta_i | \frac{\hat{\nu}_k + 1 - L}{2}, \frac{\hat{\nu}_k + 1 - L}{2})}{\sum_{i=1}^{I} \sum_{k=1}^{K} \mathbb{E}_{q_{\tau}}[\pi_k(\tau)] \mathcal{N}_D(\hat{y} | \hat{\mathbf{c}}_k^*, \hat{\mathbf{\Gamma}}_k^*(\eta_i)) \mathrm{Gam}(\eta_i | \frac{\hat{\nu}_k + 1 - L}{2}, \frac{\hat{\nu}_k + 1 - L}{2})}.
$$

Here, we used the fact that $p(\hat{y}, \hat{x} | \hat{z} = k) = p(\hat{x} | \hat{y}, \hat{z} = k)p(\hat{y} | \hat{z} = k)$, namely,

$$
\begin{align*}
& \mathcal{N}_D(\hat{y} | \hat{\mathbf{A}}_k \hat{\mathbf{x}} + \hat{\mathbf{b}}_k, \hat{\mathbf{\Sigma}}_k) \mathcal{N}_L(\hat{x} | \hat{\mathbf{m}}_k, (\eta \mathbf{L}_k)^{-1}) \\
&= \mathcal{N}_L(\hat{x} | \hat{\mathbf{\Sigma}}_k^* [\hat{\mathbf{A}}_k^\top \hat{\mathbf{\Sigma}}_k^{-1} (\hat{y} - \hat{\mathbf{b}}_k) + \eta \mathbf{L}_k \hat{\mathbf{m}}_k], \hat{\mathbf{\Sigma}}_k^*) \mathcal{N}_D(\hat{y} | \hat{\mathbf{A}}_k \hat{\mathbf{m}}_k + \hat{\mathbf{b}}_k, \hat{\mathbf{\Sigma}}_k + \hat{\mathbf{A}}_k (\eta \mathbf{L}_k)^{-1} \hat{\mathbf{A}}_k^\top) \\
&= \mathcal{N}_L(\hat{x} | \hat{\mathbf{A}}_k^*(\eta)\hat{y} + \hat{\mathbf{b}}_k^*(\eta), \hat{\mathbf{\Sigma}}_k^*(\eta)) \mathcal{N}_D(\hat{y} | \hat{\mathbf{c}}_k^*, \hat{\mathbf{\Gamma}}_k^*(\eta)),
\end{align*}
$$

where

$$
\hat{\mathbf{\Sigma}}_k^*(\eta) = (\eta \mathbf{L}_k + \hat{\mathbf{A}}_k^\top \hat{\mathbf{\Sigma}}_k^{-1} \hat{\mathbf{A}}_k)^{-1},
$$

$$
\hat{\mathbf{A}}_k^*(\eta) = \hat{\mathbf{\Sigma}}_k^*(\eta) \hat{\mathbf{A}}_k^\top \hat{\mathbf{\Sigma}}_k^{-1},
$$

$$
\hat{\mathbf{b}}_k^*(\eta) = \hat{\mathbf{\Sigma}}_k^*(\eta) [\eta \mathbf{L}_k \hat{\mathbf{m}}_k - \hat{\mathbf{A}}_k^\top \hat{\mathbf{\Sigma}}_k^{-1} \hat{\mathbf{b}}_k],
$$

$$
\hat{\mathbf{c}}_k^* = \hat{\mathbf{A}}_k \hat{\mathbf{m}}_k + \hat{\mathbf{b}}_k,
$$

$$
\hat{\mathbf{\Gamma}}_k^*(\eta) = \hat{\mathbf{\Sigma}}_k + \hat{\mathbf{A}}_k (\eta \mathbf{L}_k)^{-1} \hat{\mathbf{A}}_k^\top.
$$

The last approximation is deduced by using the fact that one simplistic strategy for evaluating integration would be to discretize $\eta$-space (1-dimensional) into a uniform grid and to evaluate the integrand as a Riemann integral with a truncated value $0 < U_{\eta} < \infty$ and a number of point $I \in \mathbb{N}^*$ for approximating the integration.
Therefore, we obtain

$p(\hat{x} | \hat{y}, \mathcal{X}, \mathcal{Y})$

$$
\begin{align*}
&\approx \sum_{k=1}^{K} \frac{\mathbb{E}_{q_{\tau}}[\pi_k(\tau)] \mathrm{St}(\hat{x} | \hat{\mathbf{m}}_k, \mathbf{L}_k, \hat{\nu}_k + 1 - L) \mathcal{N}_D(\hat{y} | \hat{\mathbf{A}}_k \hat{\mathbf{x}} + \hat{\mathbf{b}}_k, \hat{\mathbf{\Sigma}}_k)}{\sum_{k=1}^{K} \mathbb{E}_{q_{\tau}}[\pi_k(\tau)] \int \mathrm{St}(\hat{x} | \hat{\mathbf{m}}_k, \mathbf{L}_k, \hat{\nu}_k + 1 - L) \mathcal{N}_D(\hat{y} | \hat{\mathbf{A}}_k \hat{\mathbf{x}} + \hat{\mathbf{b}}_k, \hat{\mathbf{\Sigma}}_k) d\hat{\mathbf{x}}} \\
&= \sum_{k=1}^{K} \frac{\mathbb{E}_{q_{\tau}}[\pi_k(\tau)] \int_{0}^{\infty} \mathcal{N}_D(\hat{y} | \hat{\mathbf{A}}_k \hat{\mathbf{x}} + \hat{\mathbf{b}}_k, \hat{\mathbf{\Sigma}}_k) \mathcal{N}_L(\hat{x} | \hat{\mathbf{m}}_k, (\eta \mathbf{L}_k)^{-1}) \mathrm{Gam}(\eta | \frac{\hat{\nu}_k + 1 - L}{2}, \frac{\hat{\nu}_k + 1 - L}{2}) d\eta}{\sum_{k=1}^{K} \mathbb{E}_{q_{\tau}}[\pi_k(\tau)] \int_{0}^{\infty} \mathcal{N}_D(\hat{y} | \hat{\mathbf{A}}_k \hat{\mathbf{m}}_k + \hat{\mathbf{b}}_k, \hat{\mathbf{\Sigma}}_k + \eta^{-1} \hat{\mathbf{A}}_k \mathbf{L}_k^{-1} \hat{\mathbf{A}}_k^{\top}) \mathrm{Gam}(\eta | \frac{\hat{\nu}_k + 1 - L}{2}, \frac{\hat{\nu}_k + 1 - L}{2}) d\eta} \\
&= \sum_{k=1}^{K} \frac{\mathbb{E}_{q_{\tau}}[\pi_k(\tau)] \int_{0}^{\infty} \mathcal{N}_D(\hat{y} | \hat{\mathbf{c}}_k^*, \hat{\mathbf{\Gamma}}_k^*(\eta)) \mathcal{N}_L(\hat{x} | \hat{\mathbf{A}}_k^*(\eta)\hat{y} + \hat{\mathbf{b}}_k^*(\eta), \hat{\mathbf{\Sigma}}_k^*(\eta)) \mathrm{Gam}(\eta | \frac{\hat{\nu}_k + 1 - L}{2}, \frac{\hat{\nu}_k + 1 - L}{2}) d\eta}{\sum_{k=1}^{K} \mathbb{E}_{q_{\tau}}[\pi_k(\tau)] \int_{0}^{\infty} \mathcal{N}_D(\hat{y} | \hat{\mathbf{c}}_k^*, \hat{\mathbf{\Gamma}}_k^*(\eta)) \mathrm{Gam}(\eta | \frac{\hat{\nu}_k + 1 - L}{2}, \frac{\hat{\nu}_k + 1 - L}{2}) d\eta} \\
&\approx \sum_{i=1}^{I} \sum_{k=1}^{K} g_{ki}(\hat{y} | \hat{\Theta}^*, \hat{\phi}^*, \mathcal{X}, \mathcal{Y}) \mathcal{N}_L(\hat{x} | \hat{\mathbf{A}}_k^*(\eta_i)\hat{y} + \hat{\mathbf{b}}_k^*(\eta_i), \hat{\mathbf{\Sigma}}_k^*(\eta_i)),
\end{align*}
$$

where, for all $k \in [K], i \in [I]$,

$$
g_{ki}(\hat{y} | \hat{\Theta}^*, \hat{\phi}^*, \mathcal{X}, \mathcal{Y}) = \frac{\mathbb{E}_{q_{\tau}}[\pi_k(\tau)] \mathcal{N}_D(\hat{y} | \hat{\mathbf{c}}_k^*, \hat{\mathbf{\Gamma}}_k^*(\eta_i)) \mathrm{Gam}(\eta_i | \frac{\hat{\nu}_k + 1 - L}{2}, \frac{\hat{\nu}_k + 1 - L}{2})}{\sum_{i=1}^{I} \sum_{k=1}^{K} \mathbb{E}_{q_{\tau}}[\pi_k(\tau)] \mathcal{N}_D(\hat{y} | \hat{\mathbf{c}}_k^*, \hat{\mathbf{\Gamma}}_k^*(\eta_i)) \mathrm{Gam}(\eta_i | \frac{\hat{\nu}_k + 1 - L}{2}, \frac{\hat{\nu}_k + 1 - L}{2})}.
$$

Here, we used the fact that $p(\hat{y}, \hat{x} | \hat{z} = k) = p(\hat{x} | \hat{y}, \hat{z} = k)p(\hat{y} | \hat{z} = k)$, namely,

$$
\begin{align*}
& \mathcal{N}_D(\hat{y} | \hat{\mathbf{A}}_k \hat{\mathbf{x}} + \hat{\mathbf{b}}_k, \hat{\mathbf{\Sigma}}_k) \mathcal{N}_L(\hat{x} | \hat{\mathbf{m}}_k, (\eta \mathbf{L}_k)^{-1}) \\
&= \mathcal{N}_L(\hat{x} | \hat{\mathbf{\Sigma}}_k^* [\hat{\mathbf{A}}_k^\top \hat{\mathbf{\Sigma}}_k^{-1} (\hat{y} - \hat{\mathbf{b}}_k) + \eta \mathbf{L}_k \hat{\mathbf{m}}_k], \hat{\mathbf{\Sigma}}_k^*) \mathcal{N}_D(\hat{y} | \hat{\mathbf{A}}_k \hat{\mathbf{m}}_k + \hat{\mathbf{b}}_k, \hat{\mathbf{\Sigma}}_k + \hat{\mathbf{A}}_k (\eta \mathbf{L}_k)^{-1} \hat{\mathbf{A}}_k^\top) \\
&= \mathcal{N}_L(\hat{x} | \hat{\mathbf{A}}_k^*(\eta)\hat{y} + \hat{\mathbf{b}}_k^*(\eta), \hat{\mathbf{\Sigma}}_k^*(\eta)) \mathcal{N}_D(\hat{y} | \hat{\mathbf{c}}_k^*, \hat{\mathbf{\Gamma}}_k^*(\eta)),
\end{align*}
$$

where

$$
\begin{align*}
\hat{\mathbf{\Sigma}}_k^*(\eta) &= (\eta \mathbf{L}_k + \hat{\mathbf{A}}_k^\top \hat{\mathbf{\Sigma}}_k^{-1} \hat{\mathbf{A}}_k)^{-1}, \\
\hat{\mathbf{A}}_k^*(\eta) &= \hat{\mathbf{\Sigma}}_k^*(\eta) \hat{\mathbf{A}}_k^\top \hat{\mathbf{\Sigma}}_k^{-1}, \\
\hat{\mathbf{b}}_k^*(\eta) &= \hat{\mathbf{\Sigma}}_k^*(\eta) [\eta \mathbf{L}_k \hat{\mathbf{m}}_k - \hat{\mathbf{A}}_k^\top \hat{\mathbf{\Sigma}}_k^{-1} \hat{\mathbf{b}}_k], \\
\hat{\mathbf{c}}_k^* &= \hat{\mathbf{A}}_k \hat{\mathbf{m}}_k + \hat{\mathbf{b}}_k, \\
\hat{\mathbf{\Gamma}}_k^*(\eta) &= \hat{\mathbf{\Sigma}}_k + \hat{\mathbf{A}}_k (\eta \mathbf{L}_k)^{-1} \hat{\mathbf{A}}_k^\top.
\end{align*}
$$

The last approximation is deduced by using the fact that one simplistic strategy for evaluating integration would be to discretize $\eta$-space (1-dimensional) into a uniform grid and to evaluate the integrand as a Riemann integral with a truncated value $0 < U_{\eta} < \infty$ and a number of point $I \in \mathbb{N}^*$ for approximating the integration.