For the special case where \({Y}_{\mathcal {C}}|\mathcal {D}\sim \mathcal {N}({\mu }_{{Y}_{\mathcal {C}}},\text{diag}({\sigma }_{{Y}_{\mathcal {C}}}^2))\)  follows a Gaussian distribution with diagonal covariance, we can employ the development found in [1]} ([1]}) in order to obtain a closed-form solution to marginalize the output layer's uncertainty,
\(\begin{array}{ll}p(y^{(\mathtt {C})}_{\mathcal {C}}|\mathcal {D})&=\displaystyle \int p(y^{(\mathtt {C})}_{\mathcal {C}}|{y}_{\mathcal {C}})\cdot f({y}_{\mathcal {C}}|\mathcal {D})d{y}\\[-1pt]&\displaystyle =\prod _{h=1}^{\mathtt {H}}\Phi \!\left(\!(-1)^{\mathcal {C}_h}\frac{[{\mu }_{{Y}_{\mathcal {C}}}]_h}{\sqrt{\alpha ^2+[{\sigma }_{{Y}_{\mathcal {C}}}^2]_h}}\!\right),\end{array}\) 
