Adaptive sampling for Bayesian geospatial models

Abstract

Bayesian hierarchical modeling with Gaussian process random effects provides a popular approach for analyzing point-referenced spatial data. For large spatial data sets, however, generic posterior sampling is infeasible due to the extremely high computational burden in decomposing the spatial correlation matrix. In this paper, we propose an efficient algorithm—the adaptive griddy Gibbs (AGG) algorithm—to address the computational issues with large spatial data sets. The proposed algorithm dramatically reduces the computational complexity. We show theoretically that the proposed method can approximate the real posterior distribution accurately. The sufficient number of grid points for a required accuracy has also been derived. We compare the performance of AGG with that of the state-of-the-art methods in simulation studies. Finally, we apply AGG to spatially indexed data concerning building energy consumption.

Appendix: Proofs

Proof of Theorem 1

We first show that F _c (ϕ) is well defined. For θ=(σ ²,τ ²) in the likelihood function L(ϕ,θ) in Equation (3), the marginal posterior of ϕ becomes,

$$\begin{aligned} \pi(\phi\mid \boldsymbol {Y}) \propto& \pi(\phi)\int_{ \varTheta} \mathrm{L}(\phi,\boldsymbol{\theta})\pi(\boldsymbol{\theta})\, d\boldsymbol{ \theta}, \end{aligned}$$

where Θ is the support for posterior distribution of θ. Denote R(ϕ)=∫ _Θ L(ϕ,θ)π(θ) d θ. We first want to verify R(ϕ) is continuous for any ϕ∈[a,b]. Since |σ ² H(ϕ)+τ ² I _n |>τ ²ⁿ, π(θ) is product of proper priors, it then follows

• (i) The joint posterior for (ϕ,θ) is proper, so R(ϕ) is well defined.

• (ii) For any θ∈Θ and ϕ∈[a,b], L(ϕ,θ) ≤g(θ), for a function g with the property ∫ _θ∈Θ g(θ)π(θ) d θ<∞.

For any sequence ϕ _n →ϕ, continuity of L(ϕ _n ,θ) implies pointwise convergence to L(ϕ,θ), ∀θ∈Θ. This, along with (ii) also implies, by Dominated Convergence Theorem, R(ϕ _n )→R(ϕ), so R(ϕ) is continuous. For the continuous uniform prior π(ϕ), the marginal posterior cdf of ϕ turns out to be

$$\begin{aligned} F_C(\phi) =& \frac{1}{c_1}\int_{\phi_{\min}}^{\phi} R(\eta)\,d\eta,\quad \phi\in(\phi_{\min}, \phi_{\max}), \end{aligned}$$

where $c_{1} = \int_{\phi_{\min}}^{\phi_{\max}} R(\phi)\,d\phi$. F _C (ϕ)<∞, e.g., is well-defined.

For any ϵ>0, let k>(ϕ _max−ϕ _min)/ϵ and define

$$\begin{aligned} E_{j}=\biggl\{\phi: \frac{j - 1}{k}\leq F_C(\phi) < \frac{j}{k}\biggr\},\quad j = 1, \ldots, k. \end{aligned}$$

We further define $F_{D, k}(\phi) = \frac{j - 1}{k}$ for ϕ∈E _j . Apparently, for any ϕ∈(ϕ _min,ϕ _max), we have

$$0 \leq F_{C}(\phi) - F_{D, k}(\phi) \leq\frac{1}{k}. $$

Consequently, we obtain the following to complete the proof:

$$\begin{aligned} |F_{D,k}-F_C|_{TV} = & \int _{\phi_{\min}}^{\phi_{\max}} \bigl|F_{D,k}(\phi)-F_C( \phi)\bigr| \,d\phi \\ \leq& \frac{(\phi_{\max}-\phi_{\min})}{k} \\ \leq& \frac{(\phi_{\max}-\phi_{\min})}{(\phi_{\max} - \phi _{\min})/\epsilon} \\ \leq& \epsilon. \end{aligned}$$

 □