This document provides a detailed, step-by-step derivation of the full conditional distributions for all parameters in the exponential kriging model from equations (6.1) and (6.3) of Hierarchical Modeling and Analysis for Spatial Data (3rd Edition).
Special focus is given to deriving:
\[p(\boldsymbol{\beta} \mid \mathbf{Y}, \mathbf{W}, \tau^2) \propto p(\mathbf{Y} \mid \boldsymbol{\beta}, \mathbf{W}, \tau^2) \times p(\boldsymbol{\beta})\]
A hierarchical model (also called a multi-level model) organizes parameters into levels or stages. Each level depends on parameters from the level above it:
Level 0: Data (observed)
↑ depends on
Level 1: Latent Process (unobserved)
↑ depends on
Level 2: Parameters
↑ depends on
Level 3: Hyperparameters (fixed or given priors)\[\mathbf{Y} \mid \boldsymbol{\beta}, \mathbf{W}, \tau^2 \sim N(\mathbf{X}\boldsymbol{\beta} + \mathbf{W}, \tau^2 \mathbf{I})\]
\[\mathbf{W} \mid \sigma^2, \phi \sim N(\mathbf{0}, \sigma^2 H(\phi))\]
\[\boldsymbol{\beta} \sim N(\boldsymbol{\mu}_\beta, \Sigma_\beta)\]
\[\tau^2 \sim IG(a_\tau, b_\tau)\]
\[\sigma^2 \sim IG(a_\sigma, b_\sigma)\]
\[\phi \sim \text{Gamma}(a_\phi, b_\phi)\]
┌─────────────────┐
│ Hyperpriors │
│ (Fixed values) │
└────────┬────────┘
│
┌────────▼────────┐
│ Priors │
│ β, τ², σ², φ │
└────────┬────────┘
│
┌────────▼────────┐
│ Process Model │
│ W | σ², φ │
└────────┬────────┘
│
┌────────▼────────┐
│ Data Model │
│ Y | β, W, τ² │
└─────────────────┘Key Insight: Each level depends only on the level directly above it.
In a hierarchical model, the joint distribution equals the product of each variable’s conditional distribution given its parents:
\[\boxed{p(\text{all variables}) = \prod_{\text{each variable } V} p(V \mid \text{parents of } V)}\]
| Variable | Depends On (Parents) | Conditional Distribution |
|---|---|---|
| \(\mathbf{Y}\) | \(\boldsymbol{\beta}, \mathbf{W}, \tau^2\) | \(p(\mathbf{Y} \mid \boldsymbol{\beta}, \mathbf{W}, \tau^2)\) |
| \(\mathbf{W}\) | \(\sigma^2, \phi\) | \(p(\mathbf{W} \mid \sigma^2, \phi)\) |
| \(\boldsymbol{\beta}\) | None (prior) | \(p(\boldsymbol{\beta})\) |
| \(\tau^2\) | None (prior) | \(p(\tau^2)\) |
| \(\sigma^2\) | None (prior) | \(p(\sigma^2)\) |
| \(\phi\) | None (prior) | \(p(\phi)\) |
\[\boxed{p(\mathbf{Y}, \mathbf{W}, \boldsymbol{\beta}, \tau^2, \sigma^2, \phi) = p(\mathbf{Y} \mid \boldsymbol{\beta}, \mathbf{W}, \tau^2) \times p(\mathbf{W} \mid \sigma^2, \phi) \times p(\boldsymbol{\beta}) \times p(\tau^2) \times p(\sigma^2) \times p(\phi)}\]
\[p(\mathbf{Y}, \mathbf{W}, \boldsymbol{\beta}, \tau^2, \sigma^2, \phi) = p(\mathbf{Y} \mid \boldsymbol{\beta}, \mathbf{W}, \tau^2) \times p(\mathbf{W} \mid \sigma^2, \phi) \times p(\boldsymbol{\beta}) \times p(\tau^2) \times p(\sigma^2) \times p(\phi)\]
When we condition on \(\boldsymbol{\beta}\), we treat it as fixed/known. This means we divide by \(p(\boldsymbol{\beta})\):
\[p(\mathbf{Y}, \mathbf{W}, \tau^2, \sigma^2, \phi \mid \boldsymbol{\beta}) = \frac{p(\mathbf{Y}, \mathbf{W}, \boldsymbol{\beta}, \tau^2, \sigma^2, \phi)}{p(\boldsymbol{\beta})}\]
\[\boxed{p(\mathbf{Y}, \mathbf{W}, \tau^2, \sigma^2, \phi \mid \boldsymbol{\beta}) = p(\mathbf{Y} \mid \boldsymbol{\beta}, \mathbf{W}, \tau^2) \times p(\mathbf{W} \mid \sigma^2, \phi) \times p(\tau^2) \times p(\sigma^2) \times p(\phi)}\]
Variables are conditionally independent given their parents: - \(\mathbf{W}\) and \(\boldsymbol{\beta}\) are independent given \(\sigma^2, \phi\) - \(\tau^2\) and \(\boldsymbol{\beta}\) are independent (separate priors)
\[p(A, B, C) = p(A \mid B, C) \times p(B \mid C) \times p(C)\]
The probability of the sequence is the product of each step’s probability.
| Variable | Type | Depends On | Does NOT Depend On |
|---|---|---|---|
| \(\mathbf{Y}\) | Data | \(\boldsymbol{\beta}, \mathbf{W}, \tau^2\) | \(\sigma^2, \phi\) |
| \(\mathbf{W}\) | Latent | \(\sigma^2, \phi\) | \(\boldsymbol{\beta}, \tau^2\) |
| \(\boldsymbol{\beta}\) | Parameter | None (prior) | \(\mathbf{W}, \tau^2, \sigma^2, \phi\) |
| \(\tau^2\) | Parameter | None (prior) | \(\boldsymbol{\beta}, \mathbf{W}, \sigma^2, \phi\) |
| \(\sigma^2\) | Parameter | None (prior) | \(\boldsymbol{\beta}, \mathbf{W}, \tau^2, \phi\) |
| \(\phi\) | Parameter | None (prior) | \(\boldsymbol{\beta}, \mathbf{W}, \tau^2, \sigma^2\) |
\[p(\boldsymbol{\beta} \mid \mathbf{Y}, \mathbf{W}, \tau^2) \propto p(\mathbf{Y} \mid \boldsymbol{\beta}, \mathbf{W}, \tau^2) \times p(\boldsymbol{\beta})\]
For any two random variables/vectors \(A\) and \(B\):
\[P(A \mid B) = \frac{P(A \cap B)}{P(B)}\]
For continuous random variables with probability density functions:
\[p(\boldsymbol{\beta} \mid \mathbf{Y}, \mathbf{W}, \tau^2) = \frac{p(\boldsymbol{\beta}, \mathbf{Y}, \mathbf{W}, \tau^2)}{p(\mathbf{Y}, \mathbf{W}, \tau^2)}\]
Using the Chain Rule of Probability:
\[p(\boldsymbol{\beta}, \mathbf{Y}, \mathbf{W}, \tau^2) = p(\mathbf{Y} \mid \boldsymbol{\beta}, \mathbf{W}, \tau^2) \times p(\boldsymbol{\beta}, \mathbf{W}, \tau^2)\]
Now apply the Chain Rule again to \(p(\boldsymbol{\beta}, \mathbf{W}, \tau^2)\):
\[p(\boldsymbol{\beta}, \mathbf{W}, \tau^2) = p(\mathbf{W}, \tau^2 \mid \boldsymbol{\beta}) \times p(\boldsymbol{\beta})\]
Substituting back:
\[p(\boldsymbol{\beta}, \mathbf{Y}, \mathbf{W}, \tau^2) = p(\mathbf{Y} \mid \boldsymbol{\beta}, \mathbf{W}, \tau^2) \times p(\mathbf{W}, \tau^2 \mid \boldsymbol{\beta}) \times p(\boldsymbol{\beta})\]
\[p(\boldsymbol{\beta} \mid \mathbf{Y}, \mathbf{W}, \tau^2) = \frac{p(\mathbf{Y} \mid \boldsymbol{\beta}, \mathbf{W}, \tau^2) \times p(\mathbf{W}, \tau^2 \mid \boldsymbol{\beta}) \times p(\boldsymbol{\beta})}{p(\mathbf{Y}, \mathbf{W}, \tau^2)}\]
The denominator \(p(\mathbf{Y}, \mathbf{W}, \tau^2)\) is the marginal likelihood:
\[p(\mathbf{Y}, \mathbf{W}, \tau^2) = \int p(\mathbf{Y}, \mathbf{W}, \tau^2 \mid \boldsymbol{\beta}) \, p(\boldsymbol{\beta}) \, d\boldsymbol{\beta}\]
Important: This integral is over \(\boldsymbol{\beta}\). The result is a number that does not depend on \(\boldsymbol{\beta}\).
Therefore, for the purpose of finding the distribution of \(\boldsymbol{\beta}\), the denominator is a constant.
\[p(\boldsymbol{\beta} \mid \mathbf{Y}, \mathbf{W}, \tau^2) \propto p(\mathbf{Y} \mid \boldsymbol{\beta}, \mathbf{W}, \tau^2) \times p(\mathbf{W}, \tau^2 \mid \boldsymbol{\beta}) \times p(\boldsymbol{\beta})\]
In our hierarchical model:
\[\mathbf{W} \mid \sigma^2, \phi \sim N(\mathbf{0}, \sigma^2 H(\phi))\]
\[\tau^2 \sim IG(a_\tau, b_\tau)\]
Key observation: \(\mathbf{W}\) and \(\tau^2\) are independent of \(\boldsymbol{\beta}\).
Why? - \(\mathbf{W}\) depends only on \(\sigma^2\) and \(\phi\) (not on \(\boldsymbol{\beta}\)) - \(\tau^2\) has its own prior (not on \(\boldsymbol{\beta}\))
Therefore:
\[p(\mathbf{W}, \tau^2 \mid \boldsymbol{\beta}) = p(\mathbf{W} \mid \sigma^2, \phi) \times p(\tau^2)\]
Recall the full joint distribution from the hierarchical structure:
\[p(\mathbf{Y}, \mathbf{W}, \boldsymbol{\beta}, \tau^2, \sigma^2, \phi) = p(\mathbf{Y} \mid \boldsymbol{\beta}, \mathbf{W}, \tau^2) \times p(\mathbf{W} \mid \sigma^2, \phi) \times p(\boldsymbol{\beta}) \times p(\tau^2) \times p(\sigma^2) \times p(\phi)\]
Now, condition on \(\boldsymbol{\beta}\) (treat it as known):
\[p(\mathbf{Y}, \mathbf{W}, \tau^2, \sigma^2, \phi \mid \boldsymbol{\beta}) = \frac{p(\mathbf{Y}, \mathbf{W}, \boldsymbol{\beta}, \tau^2, \sigma^2, \phi)}{p(\boldsymbol{\beta})}\]
\[= p(\mathbf{Y} \mid \boldsymbol{\beta}, \mathbf{W}, \tau^2) \times p(\mathbf{W} \mid \sigma^2, \phi) \times p(\tau^2) \times p(\sigma^2) \times p(\phi)\]
Now integrate out \(\sigma^2\) and \(\phi\):
\[p(\mathbf{Y}, \mathbf{W}, \tau^2 \mid \boldsymbol{\beta}) = \int\int p(\mathbf{Y}, \mathbf{W}, \tau^2, \sigma^2, \phi \mid \boldsymbol{\beta}) \, d\sigma^2 \, d\phi\]
\[= p(\mathbf{Y} \mid \boldsymbol{\beta}, \mathbf{W}, \tau^2) \times p(\tau^2) \times \int p(\mathbf{W} \mid \sigma^2, \phi) \, p(\sigma^2) \, p(\phi) \, d\sigma^2 \, d\phi\]
The integral is just a constant \(C\) that does not depend on \(\boldsymbol{\beta}\):
\[p(\mathbf{Y}, \mathbf{W}, \tau^2 \mid \boldsymbol{\beta}) = C \times p(\mathbf{Y} \mid \boldsymbol{\beta}, \mathbf{W}, \tau^2) \times p(\tau^2)\]
\[p(\boldsymbol{\beta} \mid \mathbf{Y}, \mathbf{W}, \tau^2) \propto p(\mathbf{Y} \mid \boldsymbol{\beta}, \mathbf{W}, \tau^2) \times C \times p(\tau^2) \times p(\boldsymbol{\beta})\]
Since \(C\) and \(p(\tau^2)\) are constants with respect to \(\boldsymbol{\beta}\):
\[\boxed{p(\boldsymbol{\beta} \mid \mathbf{Y}, \mathbf{W}, \tau^2) \propto p(\mathbf{Y} \mid \boldsymbol{\beta}, \mathbf{W}, \tau^2) \times p(\boldsymbol{\beta})}\]
Step 1: Start with Bayes' Theorem
↓
p(β | Y, W, τ²) = p(Y | β, W, τ²) × p(β, W, τ²) / p(Y, W, τ²)
↓
Step 2: Apply Chain Rule
↓
p(β, W, τ²) = p(W, τ² | β) × p(β)
↓
Step 3: Substitute
↓
p(β | Y, W, τ²) = p(Y | β, W, τ²) × p(W, τ² | β) × p(β) / p(Y, W, τ²)
↓
Step 4: Identify constants
↓
p(W, τ² | β) = p(W | σ², φ) × p(τ²) [from hierarchical model]
p(Y, W, τ²) = constant (doesn't depend on β)
↓
Step 5: Drop constants
↓
p(β | Y, W, τ²) ∝ p(Y | β, W, τ²) × p(β)| Term | Meaning | Role |
|---|---|---|
| \(p(\boldsymbol{\beta} \mid \mathbf{Y}, \mathbf{W}, \tau^2)\) | Posterior | What we know about \(\boldsymbol{\beta}\) after seeing data |
| \(p(\mathbf{Y} \mid \boldsymbol{\beta}, \mathbf{W}, \tau^2)\) | Likelihood | How well \(\boldsymbol{\beta}\) explains the data \(\mathbf{Y}\) (given \(\mathbf{W}\) and \(\tau^2\)) |
| \(p(\boldsymbol{\beta})\) | Prior | What we knew about \(\boldsymbol{\beta}\) before seeing data |
The equation says:
“Our updated belief about \(\boldsymbol{\beta}\) is proportional to how well it explains the data times what we believed before seeing the data.”
If \(p(\boldsymbol{\beta}) \propto 1\) (constant):
\[p(\boldsymbol{\beta} \mid \mathbf{Y}, \mathbf{W}, \tau^2) \propto p(\mathbf{Y} \mid \boldsymbol{\beta}, \mathbf{W}, \tau^2)\]
The posterior is proportional to the likelihood — the data dominates.
If \(p(\boldsymbol{\beta})\) is very concentrated (small variance):
\[p(\boldsymbol{\beta} \mid \mathbf{Y}, \mathbf{W}, \tau^2) \approx p(\boldsymbol{\beta})\]
The prior dominates — the data has little influence.
\[Y(\mathbf{s}) = \mathbf{x}(\mathbf{s})^T \boldsymbol{\beta} + W(\mathbf{s}) + \epsilon(\mathbf{s})\]
Where: - \(W(\mathbf{s}) \sim GP(0, \sigma^2 H(\phi))\) - \(\epsilon(\mathbf{s}) \sim N(0, \tau^2)\) - \(H(\phi)_{ij} = \exp(-\phi \|\mathbf{s}_i - \mathbf{s}_j\|)\)
Stage 1:
\[\mathbf{Y} \mid \boldsymbol{\beta}, \mathbf{W}, \tau^2 \sim N(\mathbf{X}\boldsymbol{\beta} + \mathbf{W}, \tau^2 \mathbf{I})\]
Stage 2:
\[\mathbf{W} \mid \sigma^2, \phi \sim N(\mathbf{0}, \sigma^2 H(\phi))\]
\[\mathbf{Y} \mid \boldsymbol{\beta}, \sigma^2, \tau^2, \phi \sim N(\mathbf{X}\boldsymbol{\beta}, \sigma^2 H(\phi) + \tau^2 \mathbf{I})\]
\[p(\boldsymbol{\beta} \mid \mathbf{Y}, \mathbf{W}, \tau^2) \propto p(\mathbf{Y} \mid \boldsymbol{\beta}, \mathbf{W}, \tau^2) \times p(\boldsymbol{\beta})\]
| Term | Depends on \(\boldsymbol{\beta}\)? | Status |
|---|---|---|
| \(p(\mathbf{Y} \mid \boldsymbol{\beta}, \mathbf{W}, \tau^2)\) | YES | Keep |
| \(p(\mathbf{W} \mid \sigma^2, \phi)\) | NO (depends on \(\sigma^2, \phi\)) | Constant → drop |
| \(p(\tau^2)\) | NO | Constant → drop |
| \(p(\sigma^2)\) | NO | Constant → drop |
| \(p(\phi)\) | NO | Constant → drop |
| \(p(\boldsymbol{\beta})\) | YES | Keep |
From \(\mathbf{Y} \mid \boldsymbol{\beta}, \mathbf{W}, \tau^2 \sim N(\mathbf{X}\boldsymbol{\beta} + \mathbf{W}, \tau^2 \mathbf{I})\):
\[p(\mathbf{Y} \mid \boldsymbol{\beta}, \mathbf{W}, \tau^2) \propto \exp\left\{-\frac{1}{2\tau^2} (\mathbf{Y} - \mathbf{X}\boldsymbol{\beta} - \mathbf{W})^T (\mathbf{Y} - \mathbf{X}\boldsymbol{\beta} - \mathbf{W})\right\}\]
From \(\boldsymbol{\beta} \sim N(\boldsymbol{\mu}_\beta, \Sigma_\beta)\):
\[p(\boldsymbol{\beta}) \propto \exp\left\{-\frac{1}{2} (\boldsymbol{\beta} - \boldsymbol{\mu}_\beta)^T \Sigma_\beta^{-1} (\boldsymbol{\beta} - \boldsymbol{\mu}_\beta)\right\}\]
\[p(\boldsymbol{\beta} \mid \mathbf{Y}, \mathbf{W}, \tau^2) \propto \exp\left\{-\frac{1}{2\tau^2} (\mathbf{Y} - \mathbf{X}\boldsymbol{\beta} - \mathbf{W})^T (\mathbf{Y} - \mathbf{X}\boldsymbol{\beta} - \mathbf{W}) - \frac{1}{2} (\boldsymbol{\beta} - \boldsymbol{\mu}_\beta)^T \Sigma_\beta^{-1} (\boldsymbol{\beta} - \boldsymbol{\mu}_\beta)\right\}\]
This yields a multivariate normal:
\[\boxed{\boldsymbol{\beta} \mid \mathbf{Y}, \mathbf{W}, \tau^2 \sim N\left( \boldsymbol{\mu}_\beta^*, \Sigma_\beta^* \right)}\]
Where:
\[\Sigma_\beta^* = \left( \frac{1}{\tau^2} \mathbf{X}^T \mathbf{X} + \Sigma_\beta^{-1} \right)^{-1}\]
\[\boldsymbol{\mu}_\beta^* = \Sigma_\beta^* \left( \frac{1}{\tau^2} \mathbf{X}^T (\mathbf{Y} - \mathbf{W}) + \Sigma_\beta^{-1} \boldsymbol{\mu}_\beta \right)\]
If \(\Sigma_\beta^{-1} \to 0\):
\[\boldsymbol{\beta} \mid \mathbf{Y}, \mathbf{W}, \tau^2 \sim N\left( (\mathbf{X}^T \mathbf{X})^{-1} \mathbf{X}^T (\mathbf{Y} - \mathbf{W}), \tau^2 (\mathbf{X}^T \mathbf{X})^{-1} \right)\]
\[p(\mathbf{W} \mid \mathbf{Y}, \boldsymbol{\beta}, \sigma^2, \tau^2, \phi) \propto p(\mathbf{Y} \mid \boldsymbol{\beta}, \mathbf{W}, \tau^2) \times p(\mathbf{W} \mid \sigma^2, \phi)\]
\[\propto \exp\left\{-\frac{1}{2\tau^2} (\mathbf{Y} - \mathbf{X}\boldsymbol{\beta} - \mathbf{W})^T (\mathbf{Y} - \mathbf{X}\boldsymbol{\beta} - \mathbf{W}) - \frac{1}{2\sigma^2} \mathbf{W}^T H(\phi)^{-1} \mathbf{W}\right\}\]
\[\boxed{\mathbf{W} \mid \mathbf{Y}, \boldsymbol{\beta}, \sigma^2, \tau^2, \phi \sim N\left( \boldsymbol{\mu}_W^*, \Sigma_W^* \right)}\]
Where:
\[\Sigma_W^* = \left( \frac{1}{\tau^2} \mathbf{I} + \frac{1}{\sigma^2} H(\phi)^{-1} \right)^{-1}\]
\[\boldsymbol{\mu}_W^* = \Sigma_W^* \left( \frac{1}{\tau^2} (\mathbf{Y} - \mathbf{X}\boldsymbol{\beta}) \right)\]
\[\Sigma_W^* = \tau^2 \mathbf{I} - \tau^2 H(\phi) \left( \tau^2 H(\phi) + \sigma^2 \mathbf{I} \right)^{-1} \tau^2 \mathbf{I}\]
\[\boldsymbol{\mu}_W^* = H(\phi) \left( \tau^2 H(\phi) + \sigma^2 \mathbf{I} \right)^{-1} (\mathbf{Y} - \mathbf{X}\boldsymbol{\beta})\]
\[p(\tau^2 \mid \mathbf{Y}, \boldsymbol{\beta}, \mathbf{W}) \propto p(\mathbf{Y} \mid \boldsymbol{\beta}, \mathbf{W}, \tau^2) \times p(\tau^2)\]
\[\propto (\tau^2)^{-n/2} \exp\left\{-\frac{1}{2\tau^2} (\mathbf{Y} - \mathbf{X}\boldsymbol{\beta} - \mathbf{W})^T (\mathbf{Y} - \mathbf{X}\boldsymbol{\beta} - \mathbf{W})\right\} \times (\tau^2)^{-(a_\tau + 1)} \exp\left(-\frac{b_\tau}{\tau^2}\right)\]
\[\boxed{\tau^2 \mid \mathbf{Y}, \boldsymbol{\beta}, \mathbf{W} \sim IG\left( a_\tau^*, b_\tau^* \right)}\]
Where:
\[a_\tau^* = a_\tau + \frac{n}{2}\]
\[b_\tau^* = b_\tau + \frac{1}{2} (\mathbf{Y} - \mathbf{X}\boldsymbol{\beta} - \mathbf{W})^T (\mathbf{Y} - \mathbf{X}\boldsymbol{\beta} - \mathbf{W})\]
\[p(\sigma^2 \mid \mathbf{W}, \phi) \propto p(\mathbf{W} \mid \sigma^2, \phi) \times p(\sigma^2)\]
\[\propto (\sigma^2)^{-n/2} \exp\left\{-\frac{1}{2\sigma^2} \mathbf{W}^T H(\phi)^{-1} \mathbf{W}\right\} \times (\sigma^2)^{-(a_\sigma + 1)} \exp\left(-\frac{b_\sigma}{\sigma^2}\right)\]
\[\boxed{\sigma^2 \mid \mathbf{W}, \phi \sim IG\left( a_\sigma^*, b_\sigma^* \right)}\]
Where:
\[a_\sigma^* = a_\sigma + \frac{n}{2}\]
\[b_\sigma^* = b_\sigma + \frac{1}{2} \mathbf{W}^T H(\phi)^{-1} \mathbf{W}\]
\[p(\phi \mid \mathbf{W}, \sigma^2) \propto p(\mathbf{W} \mid \sigma^2, \phi) \times p(\phi)\]
\[\propto |H(\phi)|^{-1/2} \exp\left\{-\frac{1}{2\sigma^2} \mathbf{W}^T H(\phi)^{-1} \mathbf{W}\right\} \times p(\phi)\]
\(\phi\) appears inside \(H(\phi)\) in a nonlinear way through \(\exp(-\phi \|\mathbf{s}_i - \mathbf{s}_j\|)\). There is no closed-form conjugate distribution.
\[\boxed{\text{Update } \phi \text{ using Metropolis-Hastings or slice sampling}}\]
\[p(\phi \mid \mathbf{Y}, \boldsymbol{\beta}, \sigma^2, \tau^2) \propto \left| \sigma^2 H(\phi) + \tau^2 \mathbf{I} \right|^{-1/2} \exp\left\{-\frac{1}{2} (\mathbf{Y} - \mathbf{X}\boldsymbol{\beta})^T (\sigma^2 H(\phi) + \tau^2 \mathbf{I})^{-1} (\mathbf{Y} - \mathbf{X}\boldsymbol{\beta}) \right\} \times p(\phi)\]
| Parameter | Full Conditional | Distribution Type | Update Method |
|---|---|---|---|
| \(\boldsymbol{\beta}\) | \(N\left( \boldsymbol{\mu}_\beta^*, \Sigma_\beta^* \right)\) | Closed form | Gibbs |
| \(\mathbf{W}\) | \(N\left( \boldsymbol{\mu}_W^*, \Sigma_W^* \right)\) | Closed form | Gibbs |
| \(\tau^2\) | \(IG\left( a_\tau^*, b_\tau^* \right)\) | Closed form | Gibbs |
| \(\sigma^2\) | \(IG\left( a_\sigma^*, b_\sigma^* \right)\) | Closed form | Gibbs |
| \(\phi\) | Non-standard | No closed form | Metropolis-Hastings / Slice |
The marginal likelihood is obtained by integrating \(\mathbf{W}\) out of the joint distribution:
\[p(\mathbf{Y} \mid \boldsymbol{\beta}, \sigma^2, \tau^2, \phi) = \int p(\mathbf{Y} \mid \boldsymbol{\beta}, \mathbf{W}, \tau^2) \; p(\mathbf{W} \mid \sigma^2, \phi) \; d\mathbf{W}\]
Since both are normal, the integral is analytic:
\[\mathbf{Y} \mid \boldsymbol{\beta}, \sigma^2, \tau^2, \phi \sim N(\mathbf{X}\boldsymbol{\beta}, \sigma^2 H(\phi) + \tau^2 \mathbf{I})\]
After MCMC, draw \(\mathbf{W}\) from its conditional distribution:
\[p(\mathbf{W} \mid \mathbf{Y}, \boldsymbol{\beta}, \sigma^2, \tau^2, \phi) = N\left( \boldsymbol{\mu}_W^*, \Sigma_W^* \right)\]