Bayesian Tool Kit
The posterior predictive distribution is a critical concept in Bayesian inference, providing a way to predict future observations based on the observed data and the uncertainty about model parameters. It represents the distribution of a new observation $y_{\text{new}}$, given the observed data $\text{X}$, and incorporates both the likelihood of the new data and the uncertainty in the model parameters.
Definition
Mathematically, the posterior predictive distribution is given by:
$$p(y_{\text{new}} \mid \text{X}) = \int p(y_{\text{new}} \mid \theta)\, p(\theta \mid \text{X})\, d\theta$$
where:
- $p(y_{\text{new}} \mid \theta)$ is the likelihood of the new observation given the model parameters $\theta$
- $p(\theta \mid \text{X})$ is the posterior distribution of the parameters after observing the data
This distribution combines the prior information about the parameters, the likelihood of the observed data, and the uncertainty in the parameter estimates.
Purpose
The posterior predictive distribution serves several purposes:
- Model checking: It allows for the generation of simulated data from the posterior distribution, which can be compared to the observed data to assess model fit.
- Prediction of future observations: It can be used to make predictions about future data based on the current model and observed data, incorporating parameter uncertainty.
- Uncertainty quantification: It provides a full distribution for future observations rather than a point estimate, capturing uncertainty arising from both the data and the model parameters.
Derivation
We denote the observed data as $\text{X}$ and the future unobserved data as $y_{\text{new}}$. The posterior predictive distribution is derived as follows:
$$p(y_{\text{new}} \mid \text{X}) = \frac{p(y_{\text{new}},\, \text{X})}{p(\text{X})}$$
$$= \frac{\displaystyle\int p(y_{\text{new}},\, \text{X},\, \theta)\, d\theta}{p(\text{X})}$$
$$= \frac{\displaystyle\int p(y_{\text{new}} \mid \text{X}, \theta)\, p(\text{X}, \theta)\, d\theta}{p(\text{X})}$$
$$= \frac{\displaystyle\int p(y_{\text{new}} \mid \theta)\, p(\theta \mid \text{X})\, p(\text{X})\, d\theta}{p(\text{X})}$$
$$\boxed{p(y_{\text{new}} \mid \text{X}) = \int p(y_{\text{new}} \mid \theta)\, p(\theta \mid \text{X})\, d\theta}$$
The posterior predictive distribution is thus obtained by averaging the likelihood of the new observation over the posterior distribution of $\theta$, thereby fully accounting for parameter uncertainty in prediction.
Mathematical Structure of Predictive Distributions
Let us derive systematically how these distributions can be obtained using the building block models of statistical learning
- Binary outcomes with Binomial – Beta Model
- Count outcomes with Poisson – Gamma Model
- Numeric outcomes with
- Normal – Normal Model
- Normal – Gamma Model
- Normal – Normal – Gamma Model
I. BETA-BINOMIAL MODEL
1. MODEL SETUP
Observed data: $x$ successes in $n$ trials.
Likelihood:
$$x \mid \theta \sim \text{Binomial}(n,\theta)$$
$$p(x \mid \theta) = \binom{n}{x} \theta^{x} (1-\theta)^{n-x}, \quad x = 0,1,\dots,n$$
Prior:
$$\theta \sim \text{Beta}(a,b)$$
$$p(\theta) = \frac{1}{B(a,b)} \theta^{a-1}(1-\theta)^{b-1}, \quad 0 < \theta < 1$$
where
$$B(a,b) = \frac{\Gamma(a)\Gamma(b)}{\Gamma(a+b)}$$
2. PRIOR PREDICTIVE DERIVATION
Definition:
$$p(x) = \int_0^1 p(x \mid \theta)\, p(\theta)\, d\theta$$
Substitute the likelihood and the prior
$$p(x) = \int_0^1 \binom{n}{x} \theta^{x}(1-\theta)^{n-x} \cdot \frac{1}{B(a,b)} \theta^{a-1}(1-\theta)^{b-1}\, d\theta$$
Pull constants (terms not depending on $\theta$) outside the integral
$$p(x) = \frac{\binom{n}{x}}{B(a,b)} \int_0^1 \theta^{x}(1-\theta)^{n-x} \, \theta^{a-1}(1-\theta)^{b-1}\, d\theta$$
Combine powers of $\theta$ and powers of $(1-\theta)$
$$p(x) = \frac{\binom{n}{x}}{B(a,b)} \int_0^1 \theta^{x+a-1}(1-\theta)^{n-x+b-1}\, d\theta$$
Recognize the integral as the Beta function identity
$$\int_0^1 \theta^{\alpha-1}(1-\theta)^{\beta-1}\, d\theta = B(\alpha,\beta)$$ with $\alpha = x+a$ and $\beta = n-x+b$
Substitute this result
$$p(x) = \frac{\binom{n}{x}}{B(a,b)} \, B(x+a,\, n-x+b)$$
Expand using Gamma functions
$$p(x) = \binom{n}{x} \frac{\Gamma(x+a)\Gamma(n-x+b)}{\Gamma(n+a+b)} \cdot \frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)}$$
This is the PRIOR PREDICTIVE DISTRIBUTION, which is the BETA-BINOMIAL DISTRIBUTION with parameters $(n,a,b)$:
$$x \sim \text{BetaBinomial}(n,a,b)$$
$$p(x) = \binom{n}{x} \frac{B(x+a,\, n-x+b)}{B(a,b)}$$
3. POSTERIOR DERIVATION
Bayes rule:
$$p(\theta \mid x) = \frac{p(x\mid\theta)\,p(\theta)}{p(x)} \propto p(x\mid\theta)\,p(\theta)$$
Normalized posterior
$$p(\theta \mid x) = \frac{1}{B(a+x,\,b+n-x)} \theta^{a+x-1}(1-\theta)^{b+n-x-1}$$
4. POSTERIOR PREDICTIVE DERIVATION
Setup: a future experiment of $m$ new trials, with $y$ successes among them, $y = 0,1,\dots,m$.
Future likelihood:
$$p(y \mid \theta) = \binom{m}{y} \theta^{y}(1-\theta)^{m-y}$$
Definition of posterior predictive:
$$p(y \mid x) = \int_0^1 p(y \mid \theta)\, p(\theta \mid x)\, d\theta$$
Substitute the future likelihood and the posterior density
$$p(y \mid x) = \int_0^1 \binom{m}{y} \theta^{y}(1-\theta)^{m-y} \cdot \frac{1}{B(a+x,\,b+n-x)} \theta^{a+x-1}(1-\theta)^{b+n-x-1}\, d\theta$$
$$p(y \mid x) = \frac{\binom{m}{y}}{B(a+x,\,b+n-x)} \int_0^1 \theta^{y}(1-\theta)^{m-y}\,\theta^{a+x-1}(1-\theta)^{b+n-x-1}\, d\theta$$
$$p(y \mid x) = \frac{\binom{m}{y}}{B(a+x,\,b+n-x)} \int_0^1 \theta^{y+a+x-1}(1-\theta)^{m-y+b+n-x-1}\, d\theta$$
Apply the Beta function identity
$$\int_0^1 \theta^{\alpha-1}(1-\theta)^{\beta-1}\, d\theta = B(\alpha,\beta)$$ with $\alpha = y+a+x$ and $\beta = m-y+b+n-x$
Substitute this result
$$p(y \mid x) = \frac{\binom{m}{y}}{B(a+x,\,b+n-x)} \, B(y+a+x,\, m-y+b+n-x)$$
Expand using Gamma functions
$$p(y \mid x) = \binom{m}{y} \frac{\Gamma(y+a+x)\Gamma(m-y+b+n-x)}{\Gamma(m+a+x+b+n-x)} \cdot \frac{\Gamma(a+x+b+n-x)}{\Gamma(a+x)\Gamma(b+n-x)}$$
This is the POSTERIOR PREDICTIVE DISTRIBUTION, which is also a BETA-BINOMIAL DISTRIBUTION, with parameters $(m,\, a+x,\, b+n-x)$:
$$y \mid x \sim \text{BetaBinomial}(m,\, a+x,\, b+n-x)$$
$$p(y \mid x) = \binom{m}{y} \frac{B(y+a+x,\, m-y+b+n-x)}{B(a+x,\,b+n-x)}$$
II. POISSON-GAMMA MODEL
1. MODEL SETUP
Observed data: $x$ events observed over exposure (time, area, or units) $t$. For example,
- Number of disease cases observed (x) over person-years (t) of follow-up
- Number of road accidents observed over vehicle-kilometers traveled
- Number of equipment failures observed over machine operating hours
- Number of animal sightings observed over survey area (km²)
- Number of product defects observed over units produced
- Number of prediction requests served over GPU-hours
- Number of detected intrusions over network traffic volume
Likelihood:
$$x \mid \lambda \sim \text{Poisson}(\lambda t)$$
$$p(x \mid \lambda) = \frac{(\lambda t)^{x} e^{-\lambda t}}{x!}, \quad x = 0,1,2,\dots$$
Prior (rate parameterization of the Gamma):
$$\lambda \sim \text{Gamma}(a,b)$$
$$p(\lambda) = \frac{b^{a}}{\Gamma(a)} \lambda^{a-1} e^{-b\lambda}, \quad \lambda > 0$$
Gamma integral identity:
$$\int_0^\infty \lambda^{\alpha-1} e^{-\beta \lambda}\, d\lambda = \frac{\Gamma(\alpha)}{\beta^{\alpha}}$$
2. PRIOR PREDICTIVE DERIVATION
Definition:
$$p(x) = \int_0^\infty p(x \mid \lambda)\, p(\lambda)\, d\lambda$$
Substitute the likelihood and the prior
$$p(x) = \int_0^\infty \frac{(\lambda t)^{x} e^{-\lambda t}}{x!} \cdot \frac{b^{a}}{\Gamma(a)} \lambda^{a-1} e^{-b\lambda}\, d\lambda$$
$$p(x) = \frac{t^{x} b^{a}}{x!\,\Gamma(a)} \int_0^\infty \lambda^{x} e^{-\lambda t} \, \lambda^{a-1} e^{-b\lambda}\, d\lambda$$
$$p(x) = \frac{t^{x} b^{a}}{x!\,\Gamma(a)} \int_0^\infty \lambda^{x+a-1} e^{-(t+b)\lambda}\, d\lambda$$
Recognize the integral as the Gamma function identity
$$\int_0^\infty \lambda^{\alpha-1} e^{-\beta\lambda}\, d\lambda = \frac{\Gamma(\alpha)}{\beta^{\alpha}}$$ with $\alpha = x+a$ and $\beta = t+b$
$$p(x) = \frac{t^{x} b^{a}}{x!\,\Gamma(a)} \cdot \frac{\Gamma(x+a)}{(t+b)^{x+a}}$$
Rearranging,
$$p(x) = \frac{\Gamma(x+a)}{x!\,\Gamma(a)} \left(\frac{b}{t+b}\right)^{a} \left(\frac{t}{t+b}\right)^{x}$$
This is the PRIOR PREDICTIVE DISTRIBUTION, a NEGATIVE BINOMIAL DISTRIBUTION with parameters $\left(a,\, \dfrac{b}{t+b}\right)$:
$$x \sim \text{NegBinom}\left(a,\, \frac{b}{t+b}\right)$$
$$p(x) = \frac{\Gamma(x+a)}{x!\,\Gamma(a)} \left(\frac{b}{t+b}\right)^{a} \left(\frac{t}{t+b}\right)^{x}$$
3. POSTERIOR DERIVATION
Normalized posterior Refer
$$p(\lambda \mid x) = \frac{(b+t)^{a+x}}{\Gamma(a+x)} \lambda^{a+x-1} e^{-(b+t)\lambda}$$
4. POSTERIOR PREDICTIVE DERIVATION
Setup: a future observation window of exposure $m$, with $y$ events observed in it, $y = 0,1,2,\dots$, generated by the same rate $\lambda$.
Future likelihood:
$$p(y \mid \lambda) = \frac{(\lambda m)^{y} e^{-\lambda m}}{y!}$$
Definition of posterior predictive:
$$p(y \mid x) = \int_0^\infty p(y \mid \lambda)\, p(\lambda \mid x)\, d\lambda$$
Substitute the future likelihood and the posterior density
$$p(y \mid x) = \int_0^\infty \frac{(\lambda m)^{y} e^{-\lambda m}}{y!} \cdot \frac{(b+t)^{a+x}}{\Gamma(a+x)} \lambda^{a+x-1} e^{-(b+t)\lambda}\, d\lambda$$
$$p(y \mid x) = \frac{m^{y} (b+t)^{a+x}}{y!\,\Gamma(a+x)} \int_0^\infty \lambda^{y} e^{-\lambda m}\, \lambda^{a+x-1} e^{-(b+t)\lambda}\, d\lambda$$
$$p(y \mid x) = \frac{m^{y} (b+t)^{a+x}}{y!\,\Gamma(a+x)} \int_0^\infty \lambda^{y+a+x-1} e^{-(m+b+t)\lambda}\, d\lambda$$
$$\int_0^\infty \lambda^{\alpha-1} e^{-\beta\lambda}\, d\lambda = \frac{\Gamma(\alpha)}{\beta^{\alpha}}$$ with $\alpha = y+a+x$ and $\beta = m+b+t$
$$p(y \mid x) = \frac{m^{y} (b+t)^{a+x}}{y!\,\Gamma(a+x)} \cdot \frac{\Gamma(y+a+x)}{(m+b+t)^{y+a+x}}$$
$$p(y \mid x) = \frac{\Gamma(y+a+x)}{y!\,\Gamma(a+x)} \left(\frac{b+t}{m+b+t}\right)^{a+x} \left(\frac{m}{m+b+t}\right)^{y}$$
This is the POSTERIOR PREDICTIVE DISTRIBUTION, a NEGATIVE BINOMIAL DISTRIBUTION, with parameters $\left(a+x,\, \dfrac{b+t}{m+b+t}\right)$:
$$y \mid x \sim \text{NegBinom}\left(a+x,\, \frac{b+t}{m+b+t}\right)$$
$$p(y \mid x) = \frac{\Gamma(y+a+x)}{y!\,\Gamma(a+x)} \left(\frac{b+t}{m+b+t}\right)^{a+x} \left(\frac{m}{m+b+t}\right)^{y}$$
Both prior and posterior predictive distributions are Negative Binomial distributions
NOTE
If the original exposure is taken as $t=1$, the prior predictive reduces to
$$p(x) = \frac{\Gamma(x+a)}{x!\,\Gamma(a)} \left(\frac{b}{1+b}\right)^{a} \left(\frac{1}{1+b}\right)^{x}$$ and the posterior becomes $\lambda \mid x \sim \text{Gamma}(a+x,\, b+1)$, matching the commonly used single-observation Poisson-Gamma setup.
NORMAL-NORMAL MODEL WITH KNOWN VARIANCE
1. MODEL SETUP
Observed data: $n$ iid observations $x_1,\dots,x_n$, with sample mean $\bar{x}$.
Likelihood, with known variance $\sigma^2$:
$$x_i \mid \theta \sim N(\theta,\sigma^2), \quad i=1,\dots,n$$
The sample mean is a sufficient statistic:
$$\bar{x} \mid \theta \sim N\left(\theta,\, \frac{\sigma^2}{n}\right)$$
$$p(\bar{x}\mid\theta) = \frac{1}{\sqrt{2\pi \sigma^2/n}} \exp\left(-\frac{(\bar{x}-\theta)^2}{2\sigma^2/n}\right)$$
Prior:
$$\theta \sim N(\mu_0,\tau_0^2)$$
$$p(\theta) = \frac{1}{\sqrt{2\pi\tau_0^2}} \exp\left(-\frac{(\theta-\mu_0)^2}{2\tau_0^2}\right)$$
Identities used throughout:
Completing the square:
$$A\theta^2 – 2B\theta + C = A\left(\theta-\frac{B}{A}\right)^2 + \left(C-\frac{B^2}{A}\right)$$
Gaussian integral:
$$\int_{-\infty}^{\infty} \exp\left(-\frac{A}{2}(\theta-c)^2\right) d\theta = \sqrt{\frac{2\pi}{A}}$$
2. PRIOR PREDICTIVE DERIVATION
Definition:
$$p(\bar{x}) = \int_{-\infty}^{\infty} p(\bar{x}\mid\theta)\, p(\theta)\, d\theta$$
substitute the likelihood and the prior
$$p(\bar{x}) = \int_{-\infty}^{\infty} \frac{1}{\sqrt{2\pi\sigma^2/n}} \exp\left(-\frac{(\bar{x}-\theta)^2}{2\sigma^2/n}\right) \cdot \frac{1}{\sqrt{2\pi\tau_0^2}} \exp\left(-\frac{(\theta-\mu_0)^2}{2\tau_0^2}\right) d\theta$$
$$p(\bar{x}) = \frac{1}{\sqrt{2\pi\sigma^2/n}\sqrt{2\pi\tau_0^2}} \int_{-\infty}^{\infty} \exp\left(-\frac{1}{2}\left[\frac{(\bar{x}-\theta)^2}{\sigma^2/n} + \frac{(\theta-\mu_0)^2}{\tau_0^2}\right]\right) d\theta$$
expand both squared terms in $\theta$
$$\frac{(\bar{x}-\theta)^2}{\sigma^2/n} = \frac{n}{\sigma^2}\theta^2 – \frac{2n\bar{x}}{\sigma^2}\theta + \frac{n\bar{x}^2}{\sigma^2}$$
$$\frac{(\theta-\mu_0)^2}{\tau_0^2} = \frac{1}{\tau_0^2}\theta^2 – \frac{2\mu_0}{\tau_0^2}\theta + \frac{\mu_0^2}{\tau_0^2}$$
add the two expansions and collect terms in $\theta^2$, $\theta$, and constants
$$\left(\frac{n}{\sigma^2}+\frac{1}{\tau_0^2}\right)\theta^2 – 2\left(\frac{n\bar{x}}{\sigma^2}+\frac{\mu_0}{\tau_0^2}\right)\theta + \left(\frac{n\bar{x}^2}{\sigma^2}+\frac{\mu_0^2}{\tau_0^2}\right)$$
Let
$A = \frac{n}{\sigma^2}+\frac{1}{\tau_0^2}, \qquad B = \frac{n\bar{x}}{\sigma^2}+\frac{\mu_0}{\tau_0^2}, \qquad C = \frac{n\bar{x}^2}{\sigma^2}+\frac{\mu_0^2}{\tau_0^2}$ so the bracketed exponent is $A\theta^2 – 2B\theta + C$
Complete the square using the identity from Section 1
$$A\theta^2-2B\theta+C = A\left(\theta-\frac{B}{A}\right)^2 + \left(C-\frac{B^2}{A}\right)$$
substitute back into the integral, splitting the exponential
$$p(\bar{x}) = \frac{1}{\sqrt{2\pi\sigma^2/n}\sqrt{2\pi\tau_0^2}} \exp\left(-\frac{1}{2}\left(C-\frac{B^2}{A}\right)\right) \int_{-\infty}^{\infty} \exp\left(-\frac{A}{2}\left(\theta-\frac{B}{A}\right)^2\right) d\theta$$
apply the Gaussian integral identity to the remaining integral
$$\int_{-\infty}^{\infty} \exp\left(-\frac{A}{2}\left(\theta-\frac{B}{A}\right)^2\right) d\theta = \sqrt{\frac{2\pi}{A}}$$
substitute this result
$$p(\bar{x}) = \frac{1}{\sqrt{2\pi\sigma^2/n}\sqrt{2\pi\tau_0^2}} \sqrt{\frac{2\pi}{A}} \, \exp\left(-\frac{1}{2}\left(C-\frac{B^2}{A}\right)\right)$$
simplify the prefactor. Since $A = \dfrac{n\tau_0^2+\sigma^2}{\sigma^2\tau_0^2}$,
$$\sqrt{2\pi\sigma^2/n}\sqrt{2\pi\tau_0^2} = 2\pi\sqrt{\frac{\sigma^2\tau_0^2}{n}}$$
$$\frac{1}{2\pi\sqrt{\sigma^2\tau_0^2/n}} \sqrt{\frac{2\pi}{A}} = \frac{1}{\sqrt{2\pi}\sqrt{\sigma^2\tau_0^2/n}\sqrt{A}}$$
simplify $\sqrt{\sigma^2\tau_0^2/n}\cdot\sqrt{A}$
$$\sqrt{\frac{\sigma^2\tau_0^2}{n}}\cdot\sqrt{\frac{n\tau_0^2+\sigma^2}{\sigma^2\tau_0^2}} = \sqrt{\frac{n\tau_0^2+\sigma^2}{n}} = \sqrt{\tau_0^2+\frac{\sigma^2}{n}}$$
so the prefactor reduces to
$$\frac{1}{\sqrt{2\pi\left(\tau_0^2+\sigma^2/n\right)}}$$
now simplify the exponent term $C-\dfrac{B^2}{A}$.
Substituting $A = \dfrac{n\tau_0^2+\sigma^2}{\sigma^2\tau_0^2}$, $B = \dfrac{n\bar{x}\tau_0^2+\mu_0\sigma^2}{\sigma^2\tau_0^2}$, $C = \dfrac{n\bar{x}^2\tau_0^2+\mu_0^2\sigma^2}{\sigma^2\tau_0^2}$,
$$C-\frac{B^2}{A} = \frac{\left(n\bar{x}^2\tau_0^2+\mu_0^2\sigma^2\right)\left(n\tau_0^2+\sigma^2\right) – \left(n\bar{x}\tau_0^2+\mu_0\sigma^2\right)^2}{\sigma^2\tau_0^2\left(n\tau_0^2+\sigma^2\right)}$$
expand the numerator
$$\left(n\bar{x}^2\tau_0^2+\mu_0^2\sigma^2\right)\left(n\tau_0^2+\sigma^2\right) = n^2\bar{x}^2\tau_0^4 + n\bar{x}^2\tau_0^2\sigma^2 + n\mu_0^2\sigma^2\tau_0^2 + \mu_0^2\sigma^4$$
$$\left(n\bar{x}\tau_0^2+\mu_0\sigma^2\right)^2 = n^2\bar{x}^2\tau_0^4 + 2n\bar{x}\mu_0\sigma^2\tau_0^2 + \mu_0^2\sigma^4$$
subtract; the $n^2\bar{x}^2\tau_0^4$ and $\mu_0^2\sigma^4$ terms cancel
$$\text{numerator} = n\bar{x}^2\tau_0^2\sigma^2 + n\mu_0^2\sigma^2\tau_0^2 – 2n\bar{x}\mu_0\sigma^2\tau_0^2 = n\sigma^2\tau_0^2\left(\bar{x}^2+\mu_0^2-2\bar{x}\mu_0\right)$$
recognize $\bar{x}^2+\mu_0^2-2\bar{x}\mu_0 = (\bar{x}-\mu_0)^2$
$$\text{numerator} = n\sigma^2\tau_0^2(\bar{x}-\mu_0)^2$$
divide by the denominator
$$C-\frac{B^2}{A} = \frac{n\sigma^2\tau_0^2(\bar{x}-\mu_0)^2}{\sigma^2\tau_0^2\left(n\tau_0^2+\sigma^2\right)} = \frac{n(\bar{x}-\mu_0)^2}{n\tau_0^2+\sigma^2} = \frac{(\bar{x}-\mu_0)^2}{\tau_0^2+\sigma^2/n}$$
$$\therefore\quad p(\bar{x}) = \frac{1}{\sqrt{2\pi\left(\tau_0^2+\sigma^2/n\right)}} \exp\left(-\frac{(\bar{x}-\mu_0)^2}{2\left(\tau_0^2+\sigma^2/n\right)}\right)$$
This is the PRIOR PREDICTIVE DISTRIBUTION:
$$\bar{x} \sim N\left(\mu_0,\, \tau_0^2+\frac{\sigma^2}{n}\right)$$
3. POSTERIOR DERIVATION
Bayes rule:
$$p(\theta\mid\bar{x}) \propto p(\bar{x}\mid\theta)\,p(\theta)$$
substitute the likelihood and prior, dropping normalizing constants not depending on $\theta$
$$p(\theta\mid\bar{x}) \propto \exp\left(-\frac{(\bar{x}-\theta)^2}{2\sigma^2/n}\right)\exp\left(-\frac{(\theta-\mu_0)^2}{2\tau_0^2}\right)$$
combine exponents,
$$p(\theta\mid\bar{x}) \propto \exp\left(-\frac{1}{2}\left(A\theta^2-2B\theta+C\right)\right)$$ with $A,B,C$ defined exactly as in Section 2, Line 5
since $C$ does not depend on $\theta$
$$p(\theta\mid\bar{x}) \propto \exp\left(-\frac{1}{2}\left(A\theta^2-2B\theta\right)\right)$$
complete the square in $\theta$
$$A\theta^2-2B\theta = A\left(\theta-\frac{B}{A}\right)^2 – \frac{B^2}{A}$$
since $B^2/A$ does not depend on $\theta$
$$p(\theta\mid\bar{x}) \propto \exp\left(-\frac{A}{2}\left(\theta-\frac{B}{A}\right)^2\right)$$
this is the kernel of a Normal density with mean $B/A$ and variance $1/A$
$$\theta \mid \bar{x} \sim N\left(\frac{B}{A},\, \frac{1}{A}\right)$$
substitute $A = \dfrac{n}{\sigma^2}+\dfrac{1}{\tau_0^2}$ and $B = \dfrac{n\bar{x}}{\sigma^2}+\dfrac{\mu_0}{\tau_0^2}$, multiply numerator and denominator by $\sigma^2\tau_0^2$
$$\frac{B}{A} = \frac{n\bar{x}\tau_0^2+\mu_0\sigma^2}{n\tau_0^2+\sigma^2}, \qquad \frac{1}{A} = \frac{\sigma^2\tau_0^2}{n\tau_0^2+\sigma^2}$$
define the posterior parameters
$$\mu_n = \frac{n\bar{x}\tau_0^2+\mu_0\sigma^2}{n\tau_0^2+\sigma^2}, \qquad \tau_n^2 = \frac{\sigma^2\tau_0^2}{n\tau_0^2+\sigma^2}$$
The posterior:
$$\theta\mid\bar{x} \sim N(\mu_n,\tau_n^2)$$
$$p(\theta\mid\bar{x}) = \frac{1}{\sqrt{2\pi\tau_n^2}}\exp\left(-\frac{(\theta-\mu_n)^2}{2\tau_n^2}\right)$$
4. POSTERIOR PREDICTIVE DERIVATION
Setup: a future experiment of $m$ new iid observations $y_1,\dots,y_m \sim N(\theta,\sigma^2)$, with sample mean $\bar{y}$, generated by the same $\theta$.
Future likelihood:
$$\bar{y}\mid\theta \sim N\left(\theta,\frac{\sigma^2}{m}\right)$$
$$p(\bar{y}\mid\theta) = \frac{1}{\sqrt{2\pi\sigma^2/m}}\exp\left(-\frac{(\bar{y}-\theta)^2}{2\sigma^2/m}\right)$$
Definition of posterior predictive:
$$p(\bar{y}\mid\bar{x}) = \int_{-\infty}^{\infty} p(\bar{y}\mid\theta)\,p(\theta\mid\bar{x})\,d\theta$$
substitute the future likelihood and the posterior density
$$p(\bar{y}\mid\bar{x}) = \int_{-\infty}^{\infty} \frac{1}{\sqrt{2\pi\sigma^2/m}}\exp\left(-\frac{(\bar{y}-\theta)^2}{2\sigma^2/m}\right)\cdot\frac{1}{\sqrt{2\pi\tau_n^2}}\exp\left(-\frac{(\theta-\mu_n)^2}{2\tau_n^2}\right)d\theta$$
pull constants outside, combine the two exponentials
$$p(\bar{y}\mid\bar{x}) = \frac{1}{\sqrt{2\pi\sigma^2/m}\sqrt{2\pi\tau_n^2}} \int_{-\infty}^{\infty} \exp\left(-\frac{1}{2}\left[\frac{(\bar{y}-\theta)^2}{\sigma^2/m}+\frac{(\theta-\mu_n)^2}{\tau_n^2}\right]\right)d\theta$$
expand both squared terms in $\theta$
$$\frac{(\bar{y}-\theta)^2}{\sigma^2/m} = \frac{m}{\sigma^2}\theta^2 – \frac{2m\bar{y}}{\sigma^2}\theta + \frac{m\bar{y}^2}{\sigma^2}$$
$$\frac{(\theta-\mu_n)^2}{\tau_n^2} = \frac{1}{\tau_n^2}\theta^2 – \frac{2\mu_n}{\tau_n^2}\theta + \frac{\mu_n^2}{\tau_n^2}$$
add the two expansions and collect terms in $\theta^2$, $\theta$, and constants
$$\left(\frac{m}{\sigma^2}+\frac{1}{\tau_n^2}\right)\theta^2 – 2\left(\frac{m\bar{y}}{\sigma^2}+\frac{\mu_n}{\tau_n^2}\right)\theta + \left(\frac{m\bar{y}^2}{\sigma^2}+\frac{\mu_n^2}{\tau_n^2}\right)$$
define
$$A’ = \frac{m}{\sigma^2}+\frac{1}{\tau_n^2}, \qquad B’ = \frac{m\bar{y}}{\sigma^2}+\frac{\mu_n}{\tau_n^2}, \qquad C’ = \frac{m\bar{y}^2}{\sigma^2}+\frac{\mu_n^2}{\tau_n^2}$$
complete the square
$$A’\theta^2-2B’\theta+C’ = A’\left(\theta-\frac{B’}{A’}\right)^2 + \left(C’-\frac{B’^2}{A’}\right)$$
substitute back into the integral, splitting the exponential
$$p(\bar{y}\mid\bar{x}) = \frac{1}{\sqrt{2\pi\sigma^2/m}\sqrt{2\pi\tau_n^2}} \exp\left(-\frac{1}{2}\left(C’-\frac{B’^2}{A’}\right)\right) \int_{-\infty}^{\infty} \exp\left(-\frac{A’}{2}\left(\theta-\frac{B’}{A’}\right)^2\right)d\theta$$
apply the Gaussian integral identity
$$\int_{-\infty}^{\infty} \exp\left(-\frac{A’}{2}\left(\theta-\frac{B’}{A’}\right)^2\right)d\theta = \sqrt{\frac{2\pi}{A’}}$$
substitute this result
$$p(\bar{y}\mid\bar{x}) = \frac{1}{\sqrt{2\pi\sigma^2/m}\sqrt{2\pi\tau_n^2}}\sqrt{\frac{2\pi}{A’}}\,\exp\left(-\frac{1}{2}\left(C’-\frac{B’^2}{A’}\right)\right)$$
Similar to the simplification carried out in the Prior Predictive distribution (Section 2),
$$\text{prefactor: } \frac{1}{\sqrt{2\pi\left(\tau_n^2+\sigma^2/m\right)}}$$
$$\text{exponent term: } C’-\frac{B’^2}{A’} = \frac{(\bar{y}-\mu_n)^2}{\tau_n^2+\sigma^2/m}$$
substitute both results back
$$p(\bar{y}\mid\bar{x}) = \frac{1}{\sqrt{2\pi\left(\tau_n^2+\sigma^2/m\right)}} \exp\left(-\frac{(\bar{y}-\mu_n)^2}{2\left(\tau_n^2+\sigma^2/m\right)}\right)$$
This is the POSTERIOR PREDICTIVE DISTRIBUTION:
$$\bar{y}\mid\bar{x} \sim N\left(\mu_n,\, \tau_n^2+\frac{\sigma^2}{m}\right)$$
Summary
Prior predictive:
$$\bar{x} \sim N\left(\mu_0,\, \tau_0^2+\frac{\sigma^2}{n}\right)$$
parameters used: $(n,\,\mu_0,\,\tau_0^2)$
Posterior predictive:
$$\bar{y}\mid\bar{x} \sim N\left(\mu_n,\, \tau_n^2+\frac{\sigma^2}{m}\right)$$
parameters used: $(m,\,\mu_n,\,\tau_n^2)$, where
$$\mu_n = \frac{n\bar{x}\tau_0^2+\mu_0\sigma^2}{n\tau_0^2+\sigma^2}, \qquad \tau_n^2 = \frac{\sigma^2\tau_0^2}{n\tau_0^2+\sigma^2}$$
5. SPECIAL CASE
If $n=1$ and $m=1$, so a single observation $x$ is used to predict a single future observation $y$, the results reduce to
$$x \sim N\left(\mu_0,\, \tau_0^2+\sigma^2\right)$$
$$\mu_n = \frac{x\tau_0^2+\mu_0\sigma^2}{\tau_0^2+\sigma^2}, \qquad \tau_n^2 = \frac{\sigma^2\tau_0^2}{\tau_0^2+\sigma^2}$$
$$y\mid x \sim N\left(\mu_n,\, \tau_n^2+\sigma^2\right)$$
matching the commonly used single-observation Normal-Normal setup.
NORMAL-GAMMA MODEL (KNOWN MEAN)
1. MODEL SETUP
Mean $\mu$ is known. Unknown precision $\tau = 1/\sigma^2$ (precision, not variance).
Likelihood, for a single observation $x$:
$$x\mid\tau \sim N\left(\mu,\frac{1}{\tau}\right)$$
$$p(x\mid\tau) = \sqrt{\frac{\tau}{2\pi}}\,\exp\left(-\frac{\tau}{2}(x-\mu)^2\right)$$
Prior (rate parameterization of the Gamma):
$$\tau \sim \text{Gamma}(a_0,b_0)$$
$$p(\tau) = \frac{b_0^{a_0}}{\Gamma(a_0)}\tau^{a_0-1}e^{-b_0\tau}$$
Identity used:
Gamma integral:
$$\int_0^\infty \tau^{\alpha-1}e^{-\beta\tau}\,d\tau = \frac{\Gamma(\alpha)}{\beta^\alpha}$$
2. PRIOR PREDICTIVE DERIVATION
Definition:
$$p(x) = \int_0^\infty p(x\mid\tau)\,p(\tau)\,d\tau$$
Step 1: substitute the likelihood and the prior
$$p(x) = \int_0^\infty \sqrt{\frac{\tau}{2\pi}}\exp\left(-\frac{\tau}{2}(x-\mu)^2\right)\cdot\frac{b_0^{a_0}}{\Gamma(a_0)}\tau^{a_0-1}e^{-b_0\tau}\,d\tau$$
Step 2: pull constants outside the integral
$$p(x) = \frac{b_0^{a_0}}{\sqrt{2\pi}\,\Gamma(a_0)}\int_0^\infty \tau^{1/2}\,\tau^{a_0-1}\exp\left(-\frac{\tau}{2}(x-\mu)^2\right)e^{-b_0\tau}\,d\tau$$
Step 3: combine powers of $\tau$ and combine the exponentials
$$p(x) = \frac{b_0^{a_0}}{\sqrt{2\pi}\,\Gamma(a_0)}\int_0^\infty \tau^{a_0+1/2-1}\exp\left(-\tau\left[b_0+\frac{(x-\mu)^2}{2}\right]\right)d\tau$$
Step 4: apply the Gamma integral identity with $\alpha = a_0+1/2$, $\beta = b_0+\dfrac{(x-\mu)^2}{2}$
$$\int_0^\infty \tau^{\alpha-1}e^{-\beta\tau}\,d\tau = \frac{\Gamma(\alpha)}{\beta^\alpha}$$
Step 5: substitute this result
$$p(x) = \frac{b_0^{a_0}\,\Gamma(a_0+1/2)}{\sqrt{2\pi}\,\Gamma(a_0)}\left[b_0+\frac{(x-\mu)^2}{2}\right]^{-(a_0+1/2)}$$
Step 6: factor $b_0$ out of the bracket
$$b_0+\frac{(x-\mu)^2}{2} = b_0\left[1+\frac{(x-\mu)^2}{2b_0}\right]$$
Step 7: raise to the power $-(a_0+1/2)$ and combine with $b_0^{a_0}$
$$b_0^{a_0}\cdot b_0^{-(a_0+1/2)}\left[1+\frac{(x-\mu)^2}{2b_0}\right]^{-(a_0+1/2)} = b_0^{-1/2}\left[1+\frac{(x-\mu)^2}{2b_0}\right]^{-(a_0+1/2)}$$
Step 8: substitute back
$$p(x) = \frac{\Gamma(a_0+1/2)}{\sqrt{2\pi b_0}\,\Gamma(a_0)}\left[1+\frac{(x-\mu)^2}{2b_0}\right]^{-(a_0+1/2)}$$
Step 9: set $\nu = 2a_0$ and $s^2 = \dfrac{b_0}{a_0}$. Then $a_0+1/2 = (\nu+1)/2$, and the bracket exponent term rewrites as
$$\frac{(x-\mu)^2}{2b_0} = \frac{(x-\mu)^2}{\nu s^2}$$
since $\nu s^2 = 2a_0\cdot\dfrac{b_0}{a_0} = 2b_0$
Step 10: check the prefactor matches the standard Student-t normalizing constant $\dfrac{\Gamma((\nu+1)/2)}{\Gamma(\nu/2)\sqrt{\nu\pi}\,s}$
$$\frac{1}{\sqrt{\nu\pi}\,s} = \frac{1}{\sqrt{2a_0\pi}\sqrt{b_0/a_0}} = \sqrt{\frac{a_0}{2a_0\pi b_0}} = \frac{1}{\sqrt{2\pi b_0}}$$
which is exactly the prefactor obtained in Step 8
Step 11: substitute $\nu=2a_0$, $s^2=b_0/a_0$ into Step 8
$$p(x) = \frac{\Gamma((\nu+1)/2)}{\Gamma(\nu/2)\sqrt{\nu\pi}\,s}\left[1+\frac{1}{\nu}\left(\frac{x-\mu}{s}\right)^2\right]^{-(\nu+1)/2}$$
This is the density of a location-scale STUDENT-T DISTRIBUTION. The PRIOR PREDICTIVE DISTRIBUTION is:
$$x \sim t_{2a_0}\left(\mu,\ \frac{b_0}{a_0}\right)$$
with degrees of freedom $2a_0$, location $\mu$ (the known mean), and scale$^2$ $=\dfrac{b_0}{a_0}$.
3. POSTERIOR
Given $n$ observations $x_1,\dots,x_n$ with known mean $\mu$, the relevant sufficient statistic is the sum of squared deviations from the known mean:
$$S = \sum_{i=1}^n (x_i-\mu)^2$$
By conjugacy, the Gamma prior updates to:
$$\tau \mid x_1,\dots,x_n \sim \text{Gamma}(a_n,b_n)$$
$$a_n = a_0+\frac{n}{2}$$
$$b_n = b_0+\frac{1}{2}S = b_0+\frac{1}{2}\sum_{i=1}^n(x_i-\mu)^2$$
4. POSTERIOR PREDICTIVE DERIVATION
Setup: a single future observation $y$, generated by the same known mean $\mu$ and the same unknown $\tau$, predicted using the updated posterior $\text{Gamma}(a_n,b_n)$ in place of the prior.
Future likelihood:
$$y\mid\tau \sim N\left(\mu,\frac{1}{\tau}\right), \qquad p(y\mid\tau) = \sqrt{\frac{\tau}{2\pi}}\exp\left(-\frac{\tau}{2}(y-\mu)^2\right)$$
Posterior density:
$$p(\tau\mid x) = \frac{b_n^{a_n}}{\Gamma(a_n)}\tau^{a_n-1}e^{-b_n\tau}$$
Definition:
$$p(y\mid x) = \int_0^\infty p(y\mid\tau)\,p(\tau\mid x)\,d\tau$$
Step 1: substitute the future likelihood and the posterior density
$$p(y\mid x) = \int_0^\infty \sqrt{\frac{\tau}{2\pi}}\exp\left(-\frac{\tau}{2}(y-\mu)^2\right)\cdot\frac{b_n^{a_n}}{\Gamma(a_n)}\tau^{a_n-1}e^{-b_n\tau}\,d\tau$$
Step 2: pull constants outside the integral
$$p(y\mid x) = \frac{b_n^{a_n}}{\sqrt{2\pi}\,\Gamma(a_n)}\int_0^\infty \tau^{1/2}\,\tau^{a_n-1}\exp\left(-\frac{\tau}{2}(y-\mu)^2\right)e^{-b_n\tau}\,d\tau$$
Step 3: combine powers of $\tau$ and combine the exponentials
$$p(y\mid x) = \frac{b_n^{a_n}}{\sqrt{2\pi}\,\Gamma(a_n)}\int_0^\infty \tau^{a_n+1/2-1}\exp\left(-\tau\left[b_n+\frac{(y-\mu)^2}{2}\right]\right)d\tau$$
Step 4: apply the Gamma integral identity with $\alpha=a_n+1/2$, $\beta=b_n+\dfrac{(y-\mu)^2}{2}$
$$\int_0^\infty \tau^{\alpha-1}e^{-\beta\tau}\,d\tau = \frac{\Gamma(\alpha)}{\beta^\alpha}$$
Step 5: substitute this result
$$p(y\mid x) = \frac{b_n^{a_n}\,\Gamma(a_n+1/2)}{\sqrt{2\pi}\,\Gamma(a_n)}\left[b_n+\frac{(y-\mu)^2}{2}\right]^{-(a_n+1/2)}$$
Step 6: factor $b_n$ out of the bracket
$$b_n+\frac{(y-\mu)^2}{2} = b_n\left[1+\frac{(y-\mu)^2}{2b_n}\right]$$
Step 7: raise to the power $-(a_n+1/2)$ and combine with $b_n^{a_n}$
$$b_n^{a_n}\cdot b_n^{-(a_n+1/2)}\left[1+\frac{(y-\mu)^2}{2b_n}\right]^{-(a_n+1/2)} = b_n^{-1/2}\left[1+\frac{(y-\mu)^2}{2b_n}\right]^{-(a_n+1/2)}$$
Step 8: substitute back
$$p(y\mid x) = \frac{\Gamma(a_n+1/2)}{\sqrt{2\pi b_n}\,\Gamma(a_n)}\left[1+\frac{(y-\mu)^2}{2b_n}\right]^{-(a_n+1/2)}$$
Step 9: set $\nu’ = 2a_n$ and $s’^2 = \dfrac{b_n}{a_n}$. By the identical check performed in Section 2, Steps 9 and 10, this matches the standard Student-t form
$$p(y\mid x) = \frac{\Gamma((\nu’+1)/2)}{\Gamma(\nu’/2)\sqrt{\nu’\pi}\,s’}\left[1+\frac{1}{\nu’}\left(\frac{y-\mu}{s’}\right)^2\right]^{-(\nu’+1)/2}$$
This is the POSTERIOR PREDICTIVE DISTRIBUTION:
$$y\mid x \sim t_{2a_n}\left(\mu,\ \frac{b_n}{a_n}\right)$$
with degrees of freedom $2a_n$, location $\mu$ (the known mean, unchanged), and scale$^2$ $=\dfrac{b_n}{a_n}$.
5. SUMMARY
Prior predictive:
$$x \sim t_{2a_0}\left(\mu,\ \frac{b_0}{a_0}\right)$$
parameters used: $(a_0,b_0)$
Posterior predictive:
$$y\mid x \sim t_{2a_n}\left(\mu,\ \frac{b_n}{a_n}\right)$$
parameters used: $(a_n,b_n)$, from Section 3
Both are location-scale Student-t distributions with identical functional form, centered at the same known mean $\mu$ in both cases.
NORMAL-NORMAL-GAMMA MODEL
Unknown mean and unknown precision
1. MODEL SETUP
Likelihood, for a single observation $x$:
$$x \mid \mu,\tau \sim N\left(\mu,\, \frac{1}{\tau}\right)$$
$$p(x\mid\mu,\tau) = \sqrt{\frac{\tau}{2\pi}}\, \exp\left(-\frac{\tau}{2}(x-\mu)^2\right)$$
Joint Normal-Gamma prior on $(\mu,\tau)$, factored as $p(\mu,\tau) = p(\mu\mid\tau)\,p(\tau)$:
Conditional prior on $\mu$ given $\tau$:
$$\mu\mid\tau \sim N\left(\mu_0,\, \frac{1}{\kappa_0\tau}\right)$$
$$p(\mu\mid\tau) = \sqrt{\frac{\kappa_0\tau}{2\pi}}\, \exp\left(-\frac{\kappa_0\tau}{2}(\mu-\mu_0)^2\right)$$
Marginal prior on $\tau$ (rate parameterization of the Gamma):
$$\tau \sim \text{Gamma}(a_0,b_0)$$
$$p(\tau) = \frac{b_0^{a_0}}{\Gamma(a_0)}\tau^{a_0-1}e^{-b_0\tau}$$
Notation: $\mu,\tau \sim \text{NormalGamma}(\mu_0,\kappa_0,a_0,b_0)$
Identities used throughout:
Completing the square:
$$A\theta^2-2B\theta+C = A\left(\theta-\frac{B}{A}\right)^2+\left(C-\frac{B^2}{A}\right)$$
Gaussian integral:
$$\int_{-\infty}^{\infty}\exp\left(-\frac{A}{2}(\theta-c)^2\right)d\theta = \sqrt{\frac{2\pi}{A}}$$
Gamma integral:
$$\int_0^\infty \tau^{\alpha-1}e^{-\beta\tau}d\tau = \frac{\Gamma(\alpha)}{\beta^\alpha}$$
2. PRIOR PREDICTIVE DERIVATION
Definition, marginalizing over both $\mu$ and $\tau$:
$$p(x) = \int_0^\infty \int_{-\infty}^{\infty} p(x\mid\mu,\tau)\,p(\mu\mid\tau)\,p(\tau)\,d\mu\,d\tau$$
STEP A: integrate out $\mu$ first, for fixed $\tau$
Step 1: isolate the inner integral over $\mu$
$$\int_{-\infty}^{\infty} p(x\mid\mu,\tau)\,p(\mu\mid\tau)\,d\mu = \int_{-\infty}^{\infty} \sqrt{\frac{\tau}{2\pi}}\exp\left(-\frac{\tau}{2}(x-\mu)^2\right)\sqrt{\frac{\kappa_0\tau}{2\pi}}\exp\left(-\frac{\kappa_0\tau}{2}(\mu-\mu_0)^2\right)d\mu$$
Step 2: pull constants outside, combine exponentials
$$= \sqrt{\frac{\tau}{2\pi}}\sqrt{\frac{\kappa_0\tau}{2\pi}} \int_{-\infty}^{\infty} \exp\left(-\frac{\tau}{2}\left[(x-\mu)^2+\kappa_0(\mu-\mu_0)^2\right]\right)d\mu$$
Step 3: expand both squared terms in $\mu$
$$(x-\mu)^2 = \mu^2-2x\mu+x^2, \qquad \kappa_0(\mu-\mu_0)^2 = \kappa_0\mu^2-2\kappa_0\mu_0\mu+\kappa_0\mu_0^2$$
Step 4: add and collect terms in $\mu^2$, $\mu$, and constants
$$(1+\kappa_0)\mu^2 – 2(x+\kappa_0\mu_0)\mu + (x^2+\kappa_0\mu_0^2)$$
Step 5: define
$$A_\mu = 1+\kappa_0, \qquad B_\mu = x+\kappa_0\mu_0, \qquad C_\mu = x^2+\kappa_0\mu_0^2$$
Step 6: complete the square
$$A_\mu\mu^2-2B_\mu\mu+C_\mu = A_\mu\left(\mu-\frac{B_\mu}{A_\mu}\right)^2+\left(C_\mu-\frac{B_\mu^2}{A_\mu}\right)$$
Step 7: substitute back, split the exponential, apply the Gaussian integral identity to the remaining term
$$\int_{-\infty}^{\infty}\exp\left(-\frac{\tau}{2}\left[A_\mu\left(\mu-\frac{B_\mu}{A_\mu}\right)^2+\left(C_\mu-\frac{B_\mu^2}{A_\mu}\right)\right]\right)d\mu = \exp\left(-\frac{\tau}{2}\left(C_\mu-\frac{B_\mu^2}{A_\mu}\right)\right)\sqrt{\frac{2\pi}{\tau A_\mu}}$$
Step 8: simplify $C_\mu-B_\mu^2/A_\mu$
$$C_\mu-\frac{B_\mu^2}{A_\mu} = x^2+\kappa_0\mu_0^2 – \frac{(x+\kappa_0\mu_0)^2}{1+\kappa_0} = \frac{(x^2+\kappa_0\mu_0^2)(1+\kappa_0)-(x+\kappa_0\mu_0)^2}{1+\kappa_0}$$
Step 9: expand the numerator
$$(x^2+\kappa_0\mu_0^2)(1+\kappa_0) = x^2+\kappa_0x^2+\kappa_0\mu_0^2+\kappa_0^2\mu_0^2$$
$$(x+\kappa_0\mu_0)^2 = x^2+2\kappa_0x\mu_0+\kappa_0^2\mu_0^2$$
Step 10: subtract; the $x^2$ and $\kappa_0^2\mu_0^2$ terms cancel
$$\text{numerator} = \kappa_0x^2+\kappa_0\mu_0^2-2\kappa_0x\mu_0 = \kappa_0(x-\mu_0)^2$$
Step 11: so
$$C_\mu-\frac{B_\mu^2}{A_\mu} = \frac{\kappa_0(x-\mu_0)^2}{1+\kappa_0}$$
Step 12: substitute Steps 7 and 11 back into Step 2
$$\int_{-\infty}^{\infty} p(x\mid\mu,\tau)\,p(\mu\mid\tau)\,d\mu = \sqrt{\frac{\tau}{2\pi}}\sqrt{\frac{\kappa_0\tau}{2\pi}}\sqrt{\frac{2\pi}{\tau(1+\kappa_0)}}\, \exp\left(-\frac{\tau\kappa_0(x-\mu_0)^2}{2(1+\kappa_0)}\right)$$
Step 13: simplify the prefactor
$$\sqrt{\frac{\tau}{2\pi}}\sqrt{\frac{\kappa_0\tau}{2\pi}}\sqrt{\frac{2\pi}{\tau(1+\kappa_0)}} = \sqrt{\frac{\kappa_0\tau^2}{4\pi^2}}\cdot\sqrt{\frac{2\pi}{\tau(1+\kappa_0)}} = \sqrt{\frac{\kappa_0\tau}{2\pi(1+\kappa_0)}}$$
Step 14: so the $\mu$-integral collapses to
$$\int_{-\infty}^{\infty} p(x\mid\mu,\tau)\,p(\mu\mid\tau)\,d\mu = \sqrt{\frac{\kappa_0\tau}{2\pi(1+\kappa_0)}}\, \exp\left(-\frac{\kappa_0\tau(x-\mu_0)^2}{2(1+\kappa_0)}\right)$$
STEP B: integrate out $\tau$
Step 15: substitute Step 14 and $p(\tau)$ into the outer integral
$$p(x) = \int_0^\infty \sqrt{\frac{\kappa_0\tau}{2\pi(1+\kappa_0)}}\, \exp\left(-\frac{\kappa_0\tau(x-\mu_0)^2}{2(1+\kappa_0)}\right) \cdot \frac{b_0^{a_0}}{\Gamma(a_0)}\tau^{a_0-1}e^{-b_0\tau}\,d\tau$$
Step 16: pull constants outside, combine powers of $\tau$ ($\tau^{1/2}\cdot\tau^{a_0-1} = \tau^{a_0-1/2}$), combine exponentials
$$p(x) = \sqrt{\frac{\kappa_0}{2\pi(1+\kappa_0)}}\cdot\frac{b_0^{a_0}}{\Gamma(a_0)} \int_0^\infty \tau^{a_0+1/2-1}\exp\left(-\tau\left[b_0+\frac{\kappa_0(x-\mu_0)^2}{2(1+\kappa_0)}\right]\right)d\tau$$
Step 17: apply the Gamma integral identity with $\alpha = a_0+1/2$, $\beta = b_0+\dfrac{\kappa_0(x-\mu_0)^2}{2(1+\kappa_0)}$
$$\int_0^\infty \tau^{\alpha-1}e^{-\beta\tau}d\tau = \frac{\Gamma(\alpha)}{\beta^\alpha}$$
Step 18: substitute this result
$$p(x) = \sqrt{\frac{\kappa_0}{2\pi(1+\kappa_0)}}\cdot\frac{b_0^{a_0}\,\Gamma(a_0+1/2)}{\Gamma(a_0)}\cdot \left[b_0+\frac{\kappa_0(x-\mu_0)^2}{2(1+\kappa_0)}\right]^{-(a_0+1/2)}$$
Step 19: factor $b_0$ out of the bracket
$$b_0+\frac{\kappa_0(x-\mu_0)^2}{2(1+\kappa_0)} = b_0\left[1+\frac{\kappa_0(x-\mu_0)^2}{2b_0(1+\kappa_0)}\right]$$
Step 20: raise to the power $-(a_0+1/2)$ and combine with $b_0^{a_0}$
$$b_0^{a_0}\cdot b_0^{-(a_0+1/2)}\left[1+\frac{\kappa_0(x-\mu_0)^2}{2b_0(1+\kappa_0)}\right]^{-(a_0+1/2)} = b_0^{-1/2}\left[1+\frac{\kappa_0(x-\mu_0)^2}{2b_0(1+\kappa_0)}\right]^{-(a_0+1/2)}$$
Step 21: substitute back
$$p(x) = \sqrt{\frac{\kappa_0}{2\pi(1+\kappa_0)b_0}}\cdot\frac{\Gamma(a_0+1/2)}{\Gamma(a_0)}\left[1+\frac{\kappa_0(x-\mu_0)^2}{2b_0(1+\kappa_0)}\right]^{-(a_0+1/2)}$$
Step 22: set $\nu = 2a_0$ and $s^2 = \dfrac{b_0(1+\kappa_0)}{a_0\kappa_0}$. Then $a_0+1/2 = (\nu+1)/2$, and the bracket exponent term rewrites as
$$\frac{\kappa_0(x-\mu_0)^2}{2b_0(1+\kappa_0)} = \frac{(x-\mu_0)^2}{\nu s^2}$$
since $\nu s^2 = 2a_0\cdot\dfrac{b_0(1+\kappa_0)}{a_0\kappa_0} = \dfrac{2b_0(1+\kappa_0)}{\kappa_0}$, matching the denominator above after inverting.
Step 23: check the prefactor matches the standard Student-t normalizing constant $\dfrac{\Gamma((\nu+1)/2)}{\Gamma(\nu/2)\sqrt{\nu\pi}\,s}$:
$$\frac{1}{\sqrt{\nu\pi}\,s} = \frac{1}{\sqrt{2a_0\pi}\sqrt{b_0(1+\kappa_0)/(a_0\kappa_0)}} = \sqrt{\frac{a_0\kappa_0}{2a_0\pi b_0(1+\kappa_0)}} = \sqrt{\frac{\kappa_0}{2\pi b_0(1+\kappa_0)}}$$
which is exactly the prefactor obtained in Step 21.
Step 24: substitute $\nu=2a_0$, $s^2 = b_0(1+\kappa_0)/(a_0\kappa_0)$ into Step 21
$$p(x) = \frac{\Gamma((\nu+1)/2)}{\Gamma(\nu/2)\sqrt{\nu\pi}\,s}\left[1+\frac{1}{\nu}\left(\frac{x-\mu_0}{s}\right)^2\right]^{-(\nu+1)/2}$$
This is the density of a location-scale STUDENT-T DISTRIBUTION. The PRIOR PREDICTIVE DISTRIBUTION is:
$$x \sim t_{2a_0}\left(\mu_0,\ \frac{b_0(1+\kappa_0)}{a_0\kappa_0}\right)$$
with degrees of freedom $\nu=2a_0$, location $\mu_0$, and scale$^2$ $= \dfrac{b_0(1+\kappa_0)}{a_0\kappa_0}$.
3. POSTERIOR
Given $n$ observations $x_1,\dots,x_n$ with sample mean $\bar{x}$ and sum of squared deviations $s_n^2=\sum_{i=1}^n(x_i-\bar{x})^2$, the Normal-Gamma prior updates by conjugacy to:
$$\mu,\tau \mid x_1,\dots,x_n \sim \text{NormalGamma}(\mu_n,\kappa_n,a_n,b_n)$$
$$\kappa_n = \kappa_0+n$$
$$\mu_n = \frac{\kappa_0\mu_0+n\bar{x}}{\kappa_0+n}$$
$$a_n = a_0+\frac{n}{2}$$
$$b_n = b_0+\frac{1}{2}s_n^2+\frac{\kappa_0n(\bar{x}-\mu_0)^2}{2(\kappa_0+n)}$$
4. POSTERIOR PREDICTIVE DERIVATION
Setup: a single future observation $y$, generated by the same $(\mu,\tau)$, predicted using the updated posterior $\text{NormalGamma}(\mu_n,\kappa_n,a_n,b_n)$ in place of the prior.
Future likelihood:
$$y\mid\mu,\tau \sim N\left(\mu,\frac{1}{\tau}\right), \qquad p(y\mid\mu,\tau) = \sqrt{\frac{\tau}{2\pi}}\exp\left(-\frac{\tau}{2}(y-\mu)^2\right)$$
Posterior density, of identical form to Section 1 with updated parameters:
$$p(\mu\mid\tau,x) = \sqrt{\frac{\kappa_n\tau}{2\pi}}\exp\left(-\frac{\kappa_n\tau}{2}(\mu-\mu_n)^2\right), \qquad p(\tau\mid x) = \frac{b_n^{a_n}}{\Gamma(a_n)}\tau^{a_n-1}e^{-b_n\tau}$$
Definition:
$$p(y\mid x) = \int_0^\infty\int_{-\infty}^{\infty} p(y\mid\mu,\tau)\,p(\mu\mid\tau,x)\,p(\tau\mid x)\,d\mu\,d\tau$$
STEP A: integrate out $\mu$ first, for fixed $\tau$
Step 1: isolate the inner integral over $\mu$
$$\int_{-\infty}^{\infty} p(y\mid\mu,\tau)\,p(\mu\mid\tau,x)\,d\mu = \int_{-\infty}^{\infty}\sqrt{\frac{\tau}{2\pi}}\exp\left(-\frac{\tau}{2}(y-\mu)^2\right)\sqrt{\frac{\kappa_n\tau}{2\pi}}\exp\left(-\frac{\kappa_n\tau}{2}(\mu-\mu_n)^2\right)d\mu$$
Step 2: pull constants outside, combine exponentials
$$= \sqrt{\frac{\tau}{2\pi}}\sqrt{\frac{\kappa_n\tau}{2\pi}}\int_{-\infty}^{\infty}\exp\left(-\frac{\tau}{2}\left[(y-\mu)^2+\kappa_n(\mu-\mu_n)^2\right]\right)d\mu$$
Step 3: expand both squared terms in $\mu$
$$(y-\mu)^2 = \mu^2-2y\mu+y^2, \qquad \kappa_n(\mu-\mu_n)^2 = \kappa_n\mu^2-2\kappa_n\mu_n\mu+\kappa_n\mu_n^2$$
Step 4: add and collect terms in $\mu^2$, $\mu$, and constants
$$(1+\kappa_n)\mu^2-2(y+\kappa_n\mu_n)\mu+(y^2+\kappa_n\mu_n^2)$$
Step 5: define
$$A_\mu = 1+\kappa_n, \qquad B_\mu = y+\kappa_n\mu_n, \qquad C_\mu = y^2+\kappa_n\mu_n^2$$
Step 6: complete the square
$$A_\mu\mu^2-2B_\mu\mu+C_\mu = A_\mu\left(\mu-\frac{B_\mu}{A_\mu}\right)^2+\left(C_\mu-\frac{B_\mu^2}{A_\mu}\right)$$
Step 7: substitute back, split the exponential, apply the Gaussian integral identity
$$\int_{-\infty}^{\infty}\exp\left(-\frac{\tau}{2}\left[A_\mu\left(\mu-\frac{B_\mu}{A_\mu}\right)^2+\left(C_\mu-\frac{B_\mu^2}{A_\mu}\right)\right]\right)d\mu = \exp\left(-\frac{\tau}{2}\left(C_\mu-\frac{B_\mu^2}{A_\mu}\right)\right)\sqrt{\frac{2\pi}{\tau A_\mu}}$$
Step 8: simplify $C_\mu-B_\mu^2/A_\mu$, following the identical algebra of Section 2, Steps 8 to 11, with $x\to y$, $\mu_0\to\mu_n$, $\kappa_0\to\kappa_n$
$$C_\mu-\frac{B_\mu^2}{A_\mu} = \frac{\kappa_n(y-\mu_n)^2}{1+\kappa_n}$$
Step 9: substitute Steps 7 and 8 back into Step 2, and simplify the prefactor exactly as in Section 2, Steps 12 to 14
$$\int_{-\infty}^{\infty} p(y\mid\mu,\tau)\,p(\mu\mid\tau,x)\,d\mu = \sqrt{\frac{\kappa_n\tau}{2\pi(1+\kappa_n)}}\,\exp\left(-\frac{\kappa_n\tau(y-\mu_n)^2}{2(1+\kappa_n)}\right)$$
STEP B: integrate out $\tau$
Step 10: substitute Step 9 and $p(\tau\mid x)$ into the outer integral
$$p(y\mid x) = \int_0^\infty \sqrt{\frac{\kappa_n\tau}{2\pi(1+\kappa_n)}}\,\exp\left(-\frac{\kappa_n\tau(y-\mu_n)^2}{2(1+\kappa_n)}\right)\cdot\frac{b_n^{a_n}}{\Gamma(a_n)}\tau^{a_n-1}e^{-b_n\tau}\,d\tau$$
Step 11: pull constants outside, combine powers of $\tau$, combine exponentials
$$p(y\mid x) = \sqrt{\frac{\kappa_n}{2\pi(1+\kappa_n)}}\cdot\frac{b_n^{a_n}}{\Gamma(a_n)} \int_0^\infty \tau^{a_n+1/2-1}\exp\left(-\tau\left[b_n+\frac{\kappa_n(y-\mu_n)^2}{2(1+\kappa_n)}\right]\right)d\tau$$
Step 12: apply the Gamma integral identity with $\alpha=a_n+1/2$, $\beta=b_n+\dfrac{\kappa_n(y-\mu_n)^2}{2(1+\kappa_n)}$
$$p(y\mid x) = \sqrt{\frac{\kappa_n}{2\pi(1+\kappa_n)}}\cdot\frac{b_n^{a_n}\,\Gamma(a_n+1/2)}{\Gamma(a_n)}\cdot\left[b_n+\frac{\kappa_n(y-\mu_n)^2}{2(1+\kappa_n)}\right]^{-(a_n+1/2)}$$
Step 13: factor $b_n$ out of the bracket, following Section 2, Steps 19 to 21
$$p(y\mid x) = \sqrt{\frac{\kappa_n}{2\pi(1+\kappa_n)b_n}}\cdot\frac{\Gamma(a_n+1/2)}{\Gamma(a_n)}\left[1+\frac{\kappa_n(y-\mu_n)^2}{2b_n(1+\kappa_n)}\right]^{-(a_n+1/2)}$$
Step 14: set $\nu’ = 2a_n$ and $s’^2 = \dfrac{b_n(1+\kappa_n)}{a_n\kappa_n}$. By the identical check performed in Section 2, Steps 22 to 24, this matches the standard Student-t form
$$p(y\mid x) = \frac{\Gamma((\nu’+1)/2)}{\Gamma(\nu’/2)\sqrt{\nu’\pi}\,s’}\left[1+\frac{1}{\nu’}\left(\frac{y-\mu_n}{s’}\right)^2\right]^{-(\nu’+1)/2}$$
This is the POSTERIOR PREDICTIVE DISTRIBUTION:
$$y\mid x \sim t_{2a_n}\left(\mu_n,\ \frac{b_n(1+\kappa_n)}{a_n\kappa_n}\right)$$
with degrees of freedom $2a_n$, location $\mu_n$, and scale$^2$ $=\dfrac{b_n(1+\kappa_n)}{a_n\kappa_n}$.
5. SUMMARY
Prior predictive:
$$x \sim t_{2a_0}\left(\mu_0,\ \frac{b_0(1+\kappa_0)}{a_0\kappa_0}\right)$$
parameters used: $(\mu_0,\kappa_0,a_0,b_0)$
Posterior predictive:
$$y\mid x \sim t_{2a_n}\left(\mu_n,\ \frac{b_n(1+\kappa_n)}{a_n\kappa_n}\right)$$
parameters used: $(\mu_n,\kappa_n,a_n,b_n)$
Both are location-scale Student-t distributions with identical functional form.