This note relates to the paper Spatio-temporal modeling of particulate matter concentration through the SPDE approach by Cameletti et a. (Cameletti, M.; Lindgren, F.; Simpson, D.; Rue, H. Spatio-temporal modeling of particulate matter concentration through the SPDE approach. AStA Advances in Statistical Analysis 2013, 97, 109-131, doi https://doi.org/10.1007/s10182-012-0196-3.)
Here, I aim to demonstrate how the final equation for the joint posterior distribution, shown below, is derived.
\[ \begin{aligned} \pi(\theta, \xi \mid y) \propto & \; (\sigma_\varepsilon^2)^{-dT/2} \exp\left[ -\frac{1}{2\sigma_\varepsilon^2} \sum_{t=1}^T (y_t - z_t\beta - \xi_t)'(y_t - z_t\beta - \xi_t) \right] \\ & \times (\sigma_\omega^2)^{-dT/2} (1-a^2)^{d/2} |\tilde{\Sigma}|^{-T/2} \\ & \times \exp\left[ -\frac{1-a^2}{2\sigma_\omega^2} \xi_1' \tilde{\Sigma}^{-1} \xi_1 - \frac{1}{2\sigma_\omega^2} \sum_{t=2}^T (\xi_t - a\xi_{t-1})' \tilde{\Sigma}^{-1} (\xi_t - a\xi_{t-1}) \right] \\ & \times \prod_{i=1}^{\dim(\theta)} \pi(\theta_i) \end{aligned} \]
The joint posterior is:
\[ \pi(\theta, \xi \mid y) \propto \pi(y \mid \xi, \theta) \; \pi(\xi \mid \theta) \; \pi(\theta) \tag{7} \]
where \(\pi(\theta) = \prod_{i=1}^{\dim(\theta)} \pi(\theta_i)\) under independent priors.
Thus:
\[ \pi(y \mid \xi, \theta) = \prod_{t=1}^T \pi(y_t \mid \xi_t, \theta) \] \[ \pi(\xi \mid \theta) = \pi(\xi_1 \mid \theta) \prod_{t=2}^T \pi(\xi_t \mid \xi_{t-1}, \theta) \]
Substitute into (7):
\[ \pi(\theta, \xi \mid y) \propto \left[ \prod_{t=1}^T \pi(y_t \mid \xi_t, \theta) \right] \, \pi(\xi_1 \mid \theta) \, \prod_{t=2}^T \pi(\xi_t \mid \xi_{t-1}, \theta) \, \pi(\theta) \tag{8} \]
The observation equation (5):
\[ y_t = z_t \beta + \xi_t + \varepsilon_t, \quad \varepsilon_t \stackrel{\text{i.i.d.}}{\sim} N(0, \sigma_\varepsilon^2 I_d) \]
Hence:
\[ \pi(y_t \mid \xi_t, \theta) \propto (\sigma_\varepsilon^2)^{-d/2} \exp\left[ -\frac{1}{2\sigma_\varepsilon^2} (y_t - z_t\beta - \xi_t)'(y_t - z_t\beta - \xi_t) \right] \]
The state equation (6):
\[ \xi_t = a \xi_{t-1} + \omega_t, \quad \omega_t \stackrel{\text{i.i.d.}}{\sim} N(0, \sigma_\omega^2 \tilde{\Sigma}) \]
For \(t \ge 2\):
\[ \pi(\xi_t \mid \xi_{t-1}, \theta) \propto (\sigma_\omega^2)^{-d/2} |\tilde{\Sigma}|^{-1/2} \exp\left[ -\frac{1}{2\sigma_\omega^2} (\xi_t - a\xi_{t-1})' \tilde{\Sigma}^{-1} (\xi_t - a\xi_{t-1}) \right] \]
Initial condition (\(t=1\)):
\[ \xi_1 \sim N\left(0, \frac{\sigma_\omega^2}{1-a^2} \tilde{\Sigma}\right) \]
Thus:
\[ \pi(\xi_1 \mid \theta) \propto (\sigma_\omega^2)^{-d/2} (1-a^2)^{d/2} |\tilde{\Sigma}|^{-1/2} \exp\left[ -\frac{1-a^2}{2\sigma_\omega^2} \xi_1' \tilde{\Sigma}^{-1} \xi_1 \right] \]
\[ \prod_{t=1}^T \pi(y_t \mid \xi_t, \theta) \propto (\sigma_\varepsilon^2)^{-dT/2} \exp\left[ -\frac{1}{2\sigma_\varepsilon^2} \sum_{t=1}^T (y_t - z_t\beta - \xi_t)'(y_t - z_t\beta - \xi_t) \right] \]
\[ \pi(\xi_1 \mid \theta) \propto (\sigma_\omega^2)^{-d/2} (1-a^2)^{d/2} |\tilde{\Sigma}|^{-1/2} \exp\left[ -\frac{1-a^2}{2\sigma_\omega^2} \xi_1' \tilde{\Sigma}^{-1} \xi_1 \right] \]
\[ \prod_{t=2}^T \pi(\xi_t \mid \xi_{t-1}, \theta) \propto (\sigma_\omega^2)^{-d(T-1)/2} |\tilde{\Sigma}|^{-(T-1)/2} \exp\left[ -\frac{1}{2\sigma_\omega^2} \sum_{t=2}^T (\xi_t - a\xi_{t-1})' \tilde{\Sigma}^{-1} (\xi_t - a\xi_{t-1}) \right] \]
Gather powers of \(\sigma_\varepsilon^2\):
\[ (\sigma_\varepsilon^2)^{-dT/2} \]
Gather powers of \(\sigma_\omega^2\):
\[ (\sigma_\omega^2)^{-d/2} \cdot (\sigma_\omega^2)^{-d(T-1)/2} = (\sigma_\omega^2)^{-dT/2} \]
Gather powers of \(|\tilde{\Sigma}|\):
\[ |\tilde{\Sigma}|^{-1/2} \cdot |\tilde{\Sigma}|^{-(T-1)/2} = |\tilde{\Sigma}|^{-T/2} \]
Include \((1-a^2)^{d/2}\) from the initial density.
Collect exponential terms:
\[ \exp\left[ -\frac{1}{2\sigma_\varepsilon^2} \sum_{t=1}^T (y_t - z_t\beta - \xi_t)'(y_t - z_t\beta - \xi_t) - \frac{1-a^2}{2\sigma_\omega^2} \xi_1' \tilde{\Sigma}^{-1} \xi_1 - \frac{1}{2\sigma_\omega^2} \sum_{t=2}^T (\xi_t - a\xi_{t-1})' \tilde{\Sigma}^{-1} (\xi_t - a\xi_{t-1}) \right] \]
Finally, multiply by the prior product:
\[ \prod_{i=1}^{\dim(\theta)} \pi(\theta_i) \]
\[ \begin{aligned} \pi(\theta, \xi \mid y) \propto & \; (\sigma_\varepsilon^2)^{-dT/2} \exp\left[ -\frac{1}{2\sigma_\varepsilon^2} \sum_{t=1}^T (y_t - z_t\beta - \xi_t)'(y_t - z_t\beta - \xi_t) \right] \\ & \times (\sigma_\omega^2)^{-dT/2} (1-a^2)^{d/2} |\tilde{\Sigma}|^{-T/2} \\ & \times \exp\left[ -\frac{1-a^2}{2\sigma_\omega^2} \xi_1' \tilde{\Sigma}^{-1} \xi_1 - \frac{1}{2\sigma_\omega^2} \sum_{t=2}^T (\xi_t - a\xi_{t-1})' \tilde{\Sigma}^{-1} (\xi_t - a\xi_{t-1}) \right] \\ & \times \prod_{i=1}^{\dim(\theta)} \pi(\theta_i) \end{aligned} \]
The authors wrote:
“We now create a SPDE model object for a Mat´ern-like spatial covariance function using the function inla.spde2.matern(.) specifying the obtained triangulation (given by the mesh object) and the parameter \(\alpha\) = 2 and, as noted at the end of Section 3.2, it follows that the smoothness parameter \(\nu\) of the Mat´ern covariance function is equal to 1.”
The parameter \(\alpha\) was set to 2 because of the following mathematical relationship:
In the INLA framework, the parameter alpha
is used to define the Matérn covariance model through the stochastic
partial differential equation (SPDE) approach, where:
\[ (\kappa^2 - \Delta)^{\alpha/2} (\tau u) = \mathcal{W} \]
Here, \(\mathcal{W}\) is Gaussian white noise, \(\Delta\) is the Laplacian, and \(\kappa\) and \(\tau\) are scaling parameters. The relationship between \(\alpha\) and the smoothness parameter \(\nu\) of the Matérn covariance function is:
\[ \nu = \alpha - \frac{d}{2} \]
where \(d\) is the spatial dimension.
1. Dimension \(d =
2\)
We are working with a spatial mesh in 2D (given by the mesh
object in INLA, typically for a spatial domain in \(\mathbb{R}^2\)).
2. The smoothness parameter \(\nu\)
The text explicitly states:
“the smoothness parameter \(\nu\)
of the Matérn covariance function is equal to 1”
Plug these values into the relationship:
\[ \begin{aligned} \nu &= \alpha - \frac{d}{2} \\ 1 &= \alpha - \frac{2}{2} \\ 1 &= \alpha - 1 \\ \alpha &= 2 \end{aligned} \]
\(\alpha\) = 2 ensures \(\nu\) = 1 in 2D spatial models.
In many spatial statistics applications, \(\nu\) = 1 is a common choice because: - It corresponds to a mean-square differentiable process - It provides sufficient smoothness for many environmental and ecological applications - It is computationally convenient in the SPDE/INLA setting because \(\alpha\) = 2 leads to a simpler second-order differential operator rather than fractional powers