If \((x,y)\) has a multivariate Normal-Gamma distribution with parameters \(\mu\), \(P\),
\(\alpha\), and \(\beta\), then the marginal distribution of \(y\) has a Gamma
distribution with parameters \(\alpha\), \(\beta\), and conditionally on \(y\),
\(x\) has a multivariate normal distribution with expectation \(\mu\) and
precision matrix \(yP\). The probability density is proportional to
$$
f(x,y)=y^{\alpha+k/2-1}\exp(-y(\beta + (x-\mu)^tP(x-\mu)/2))
$$
where \(k\) is the dimension.