is satisfied by Fubini’s theorem.

Parameter Estimation
--------------------

With the above prior specification, the integrated posterior of ${\boldsymbol \phi}$ given ${\mathbf{y}}$ is given by $$\begin{aligned}
 \label{eqn: posterior}
    \pi({\boldsymbol \phi}\mid {\mathbf{y}}) \propto  \prod_{t=1}^s |{{\mathbf{R}}}_t|^{-1/2} |{{\mathbf{X}}}_t^\top {{\mathbf{R}}}_t^{-1} {{\mathbf{X}}}_t|^{-1/2} \{S^2({\boldsymbol \phi}_t)\}^{-(({n}_t - q_t)/2 + a_t-1)} \pi({\boldsymbol \phi}_t),\end{aligned}$$ where $\pi({\boldsymbol \phi}_t)$ refers to independent reference prior and independent Jeffreys priors. Inference based on this posterior distribution can be made in a fully Bayesian paradigm via Markov chain Monte Carlo methods. Although uncertainties in all parameters can be taken into consideration in a fully Bayesian approach, the associated computation can be too expensive in practice due to repeated evaluation of the integrated likelihood .

In what follows, we focus on empirical Bayesian inference by maximizing the posterior  to obtain the estimate of ${\boldsymbol \phi}$. In fact, the numerical optimization can be performed for each ${\boldsymbol \phi}_t$, that is, for $t=1, \ldots, s$, $$\begin{aligned}
    \hat{{\boldsymbol \phi}}_t := \underset{{\boldsymbol \phi}_t}{\text{argmax}} \left\{-\frac{1}{2}\ln |{\mathbf{R}}_t| - \frac{1}{2}\ln|{\mathbf{X}}_t^\top {\mathbf{R}}_t^{-1} {\mathbf{X}}_t| - \left(\frac{n_t-q_t}{2}+a_t-1\right)\ln S^2({\boldsymbol \phi}_t) + \ln \pi({\boldsymbol \phi}_t)\right\},\end{aligned}$$ where the maximization step can be performed using standard optimization algorithms such as the Nelder-Mead algorithm [@Nocedal2006]. Once $\hat{{\boldsymbol \phi}}_t$ is obtained, the cokriging predictions and cokriging variances can be obtained based on 