Bayesian Spatial Field Reconstruction with Unknown Distortions in Sensor Networks

TL;DR

This paper introduces robust Bayesian methods for spatial field reconstruction in sensor networks with unknown distortions, effectively handling biased and miscalibrated sensor data.

Contribution

It develops a mixture of Gaussian process experts model and two Bayesian algorithms, S-BLUE and empirical Bayes, for accurate spatial field estimation with distorted sensor data.

Findings

The methods outperform traditional estimators in distorted sensor scenarios.

Distributed algorithms significantly improve computational efficiency.

Simulation results validate robustness and accuracy of proposed approaches.

Abstract

Spatial regression of random fields based on potentially biased sensing information is proposed in this paper. One major concern in such applications is that since it is not known a-priori what the accuracy of the collected data from each sensor is, the performance can be negatively affected if the collected information is not fused appropriately. For example, the data collector may measure the phenomenon inappropriately, or alternatively, the sensors could be out of calibration, thus introducing random gain and bias to the measurement process. Such readings would be systematically distorted, leading to incorrect estimation of the spatial field. To combat this detrimental effect, we develop a robust version of the spatial field model based on a mixture of Gaussian process experts. We then develop two different approaches for Bayesian spatial field reconstruction: the first algorithm is the Spatial Best Linear Unbiased Estimator (S-BLUE), in which one considers the quadratic loss function and restricts the estimator to the linear family of transformations; the second algorithm is based on empirical Bayes, which utilises a two-stage estimation procedure to produce accurate predictive inference in the presence of "misbehaving" sensors. In addition, we develop the distributed version of these two approaches to drastically improve the computational efficiency in large-scale settings. We present extensive simulation results using both synthetic datasets and semi-synthetic datasets with real temperature measurements and simulated distortions to draw useful conclusions regarding the performance of each of the algorithms.

Tables

Table 1. TABLE I : Synthetic experiment 1 & 2 – maximum absolute deviation from the average relative MSE of the five methods.

Method	Experiment 1	Experiment 2
oracle	0.0278	0.1615
naive	0.1566	0.2064
S-BLUE	0.0410	0.1800
ICM	0.0937	0.1935
CEM	0.0908	0.1928

Table 2. TABLE II : Summary statistics of the temperature measurements.

mean	variance	skewness	kurtosis
19.7930	119.9365	-0.5429	2.9377
min	1st-quartile	median	3rd-quartile	max
-31.6667	12.7778	21.1111	28.3333	55.0000

Table 3. TABLE III : Estimated GP hyperparameters

signal mean (°C)	signal std. dev. (°C)
30.5034	4.6587
noise std. dev. (°C)	length-scale (km)
1.6340	174.3699

Table 4. TABLE IV : Experiment with real data – maximum absolute deviation from the average relative MSE of the eight methods.

Method	Centralized	Distributed
oracle	0.0728	-
naive	0.2015	-
S-BLUE	0.1132	0.1128
ICM	0.2676	0.3188
CEM	0.1838	0.1104

Equations

$$\begin{array}{rl}& \mathcal{F}:=\{f\left(\cdot \right):\mathcal{X}\mapsto \mathbb{R}\text{s.t.}f\left(\cdot \right)\sim \mathcal{G}\mathcal{P}\left(\mu \left(\cdot \right),\mathcal{C}\left(\cdot ,\cdot \right)\right),\\ & \mathrm{with}\mu \left(\mathbf{x}\right):=\mathbb{E}\left[f\left(\mathbf{x}\right)\right]:\mathcal{X}\mapsto \mathbb{R},\\ & \mathcal{C}\left(\mathbf{x},{\mathbf{x}}^{\prime }\right):=\mathbb{E}\left[\left(f\left(\mathbf{x}\right)-\mu \left(\mathbf{x}\right)\right)\left(f\left({\mathbf{x}}^{\prime }\right)-\mu \left({\mathbf{x}}^{\prime }\right)\right)\right]:\mathcal{X}\times \mathcal{X}\mapsto \mathbb{R}\},\end{array}$$

$$\begin{array}{rl}{\tilde{y}}_{n,m}& =f\left({\mathbf{x}}_n\right)+{ϵ}_{n,m}(\text{acquisition})\\ y_{n,m}& =\mathcal{T}\left({\tilde{y}}_{n,m};{\mathbf{\psi }}_n\right)(\text{distortion})\end{array}$$

$$\mathcal{T}\left(u;{\mathbf{\psi }}_n={\left(a_n,b_n\right)}^T\right):=a_{n}u+b_n,$$

$$F\sim \mathcal{G}\mathcal{P}(\mu (\cdot ),\mathcal{C}(\cdot ,\cdot )).$$

$${ϵ}_{n,m}\stackrel{\text{i.i.d.}}{\sim }\mathcal{N}(0,{\varsigma }^2).$$

$$\pi (Z_n=k)=q_k^{(n)},k=0,\dots ,K,$$

$$\left({\mathbf{\psi }}_n|Z_n=k\right)\sim {\pi }_k,$$

$${\pi }_0({\mathbf{\psi }}_n)={\delta }_{{\mathbf{\psi }}^0}.$$

$$\begin{array}{rl}\pi \left({\mathbf{\psi }}_n\right)=& q_0^{(n)}{\delta }_{{\mathbf{\psi }}^0}+\sum _{k=1}^K{q}_k^{(n)}{\pi }_k({\mathbf{\psi }}_n).\end{array}$$

$$\pi (\mathbf{\psi })=\prod _{n=1}^N\pi \left({\mathbf{\psi }}_n\right).$$

$$\begin{array}{rl}& \log p(\mathbf{\psi }|\mathbf{y})\\ =& -\frac{1}{2}[\mathrm{tr}(\mathbf{M})\log 2\pi +\mathrm{tr}(\mathbf{M}\log ({\varsigma }^2{\mathbf{A}}^2))-\log |{\varsigma }^2{\mathbf{M}}^{-1}|\\ & +\log |\mathbf{{\rm Y}}|+{\varsigma }^{-2}{1}^T{\mathbf{A}}^{-2}\mathbf{s}-{\varsigma }^{-2}{\mathbf{g}}^T{\mathbf{M}}^{-1}{\mathbf{A}}^{-2}\mathbf{g}\\ & +(\tilde{\mathbf{g}}-\mathbf{\mu }{)}^T{\mathbf{{\rm Y}}}^{-1}(\tilde{\mathbf{g}}-\mathbf{\mu })]+\log \pi (\mathbf{\psi })-\log p(\mathbf{y}).\end{array}$$

$$p(f_{\ast }|\mathbf{y})=\int p(f_{\ast }|\mathbf{y},\mathbf{\psi })p(\mathbf{\psi }|\mathbf{y})d\mathbf{\psi }.$$

$$p\left(f_{\ast }|\mathbf{y},\mathbf{\psi }\right)=$$

$${\bar{f}}_{\ast }=$$

$${\sigma }_{\ast }^2=$$

$$R[\Pi ,h]=\mathbb{E}[l(h(\mathbf{y}),f_{\ast })].$$

$$\begin{array}{rl}{\hat{h}}_{\text{S-BLUE}}=& \underset{h\in \mathcal{H}}{\mathrm{argmin}}R[\Pi ,h]=\underset{h\in \mathcal{H}}{\mathrm{argmin}}\mathbb{E}\left[{\left(h(\mathbf{y})-f_{\ast }\right)}^2\right],\end{array}$$

$$\begin{array}{rl}{\hat{h}}_{\text{S-BLUE}}(\mathbf{y})=& {\mu }_{\ast }+\mathrm{Cov}[\bar{\mathbf{g}},f_{\ast }{]}^T\mathrm{Cov}[\bar{\mathbf{g}}{]}^{-1}(\bar{\mathbf{g}}-\mathbb{E}[\bar{\mathbf{g}}]),\end{array}$$

$$\mathbb{E}[\bar{\mathbf{g}}]=$$

$$\mathrm{Cov}[\bar{\mathbf{g}},f_{\ast }]=$$

$$\begin{array}{rl}\mathrm{Cov}[\bar{\mathbf{g}}]=& \mathbb{E}[\mathbf{a}{\mathbf{a}}^T]\odot \left(\mathcal{C}+{\varsigma }^2{\mathbf{M}}^{-1}+\mathbf{\mu }{\mathbf{\mu }}^T\right)\\ & +\mathrm{diag}(\mathbf{\mu })\left(\mathbb{E}[\mathbf{a}{\mathbf{b}}^T]+\mathbb{E}[\mathbf{a}{\mathbf{b}}^T{]}^T\right)\\ & +\mathbb{E}[\mathbf{b}{\mathbf{b}}^T]-\mathbb{E}[\bar{\mathbf{g}}]\mathbb{E}[\bar{\mathbf{g}}{]}^T.\end{array}$$

$$\begin{array}{rl}R[\Pi ,{\hat{h}}_{\text{S-BLUE}}]=& {\mathcal{C}}_{\ast }-\mathrm{Cov}[\bar{\mathbf{g}},f_{\ast }{]}^T\mathrm{Cov}[\bar{\mathbf{g}}{]}^{-1}\mathrm{Cov}[\bar{\mathbf{g}},f_{\ast }].\end{array}$$

$${\hat{h}}_{\text{EB}}(\mathbf{y},\hat{\mathbf{\psi }})=\underset{h(\mathbf{y})}{\mathrm{argmin}}\mathbb{E}[l(h(\mathbf{y}),f_{\ast })|\mathbf{y},\hat{\mathbf{\psi }}].$$

$$\begin{array}{rl}{\hat{h}}_{\text{MMSE}}(\mathbf{y})=& \mathbb{E}[f_{\ast }|\mathbf{y}]={\int }_{\mathbb{R}}{f}_{\ast }p(f_{\ast }|\mathbf{y})d{f}_{\ast },\end{array}$$

$$\begin{array}{rl}{\hat{h}}_{\text{EB-MMSE}}(\mathbf{y},\hat{\mathbf{\psi }})=& \mathbb{E}[f_{\ast }|\mathbf{y},\hat{\mathbf{\psi }}]={\int }_{\mathbb{R}}{f}_{\ast }p(f_{\ast }|\mathbf{y},\hat{\mathbf{\psi }})d{f}_{\ast }\\ \approx & {\hat{h}}_{\text{MMSE}}(\mathbf{y}),\end{array}$$

$$\begin{array}{rl}{\hat{h}}_{\text{LAD}}(\mathbf{y})=& \mathrm{median}(f_{\ast }|\mathbf{y}).\end{array}$$

$$\begin{array}{rl}{\hat{h}}_{\text{EB-LAD}}(\mathbf{y},\hat{\mathbf{\psi }})=& \mathrm{median}(f_{\ast }|\mathbf{y},\hat{\mathbf{\psi }})\approx {\hat{h}}_{\text{LAD}}(\mathbf{y}).\end{array}$$

$$\begin{array}{rl}{\hat{h}}_{\text{MAP}}(\mathbf{y})=& \underset{f_{\ast }}{\mathrm{argmax}}p(f_{\ast }|\mathbf{y}).\end{array}$$

$${\hat{h}}_{\text{EB-MAP}}(\mathbf{y},\hat{\mathbf{\psi }})=$$

$$\displaystyle\begin{split}\hat{{\bm{\psi}}}=&\operatorname*{arg\,max}_{{\bm{\psi}}}p\left({\bm{\psi}}|{\bm{y}}\right)\\ =&\operatorname*{arg\,max}_{{\bm{\psi}}}\left[-\frac{1}{2}\left\{\mathrm{tr}(\mathbf{M})\log 2\pi+\mathrm{tr}(\mathbf{M}\log(\varsigma^{2}\mathbf{A}^{2}))\right.\right.\\ &\left.-\log|\varsigma^{2}\mathbf{M}^{-1}|+\log|{\bm{\Upsilon}}|+\varsigma^{-2}\mathbf{1}^{T}\mathbf{A}^{-2}{\bm{s}}\right.\\ &\left.-\varsigma^{-2}{\bm{g}}^{T}\mathbf{M}^{-1}\mathbf{A}^{-2}{\bm{g}}+(\tilde{{\bm{g}}}-{\bm{\mu}})^{T}{\bm{\Upsilon}}^{-1}(\tilde{{\bm{g}}}-{\bm{\mu}})\right\}\\ &+\log\pi({\bm{\psi}})\bigg{]}.\end{split}$$

$$\hat{\mathbf{\phi }}=\underset{\mathbf{\phi }\in \Phi }{\mathrm{argmax}}J(\mathbf{\phi }),$$

$${\mathbb{P}}_{\mathbf{\theta }}[J(\mathbf{\phi })\ge \gamma ]={\mathbb{E}}_{\mathbf{\theta }}[\text{mathbbm}{1}_{\{J(\mathbf{\phi })\ge \gamma \}}]={\int }_{\Phi }\text{mathbbm}{1}_{\{J(\mathbf{\phi })\ge \gamma \}}{w}_{\mathbf{\theta }}(\mathbf{\phi })d\mathbf{\phi }.$$

Full text

**Keywords: ** Sensor Networks, Gaussian Process, Spatial Linear Unbiased Estimator (S-BLUE), Empirical Bayes, Cross Entropy method (CEM), Iterated Conditional Modes (ICM)

I Introduction

In recent years, Wireless Sensor Networks (WSNs) have attracted considerable attention due to their applications in environment monitoring [1, 2], forecasting [3], surveillance [4], event detection [5] and tracking [6]. For example, the United States Environmental Protection Agency (EPA) proposed to promote the use of sensor networks for air quality monitoring [7]. In this paper, we focus on environmental monitoring applications in which a WSN consists of a collection of spatially distributed sensor nodes with limited energy and communication bandwidth. The sensors make observations of spatial physical phenomena (e.g. the concentration of air pollutants such as carbon monoxide and ozone, temperature, humidity, etc. [8, 2]) and communicate the observations to a Fusion Center (FC) [9]. The FC then reconstructs the spatial phenomena from these observations at any spatial location of interest, based on which decisions can be made and actions can be performed.

In many cases, the sensor nodes used to collect the spatial information are unreliable and distort the spatial information. These distortions should be accounted for in spatial field reconstruction. It is therefore crucial to assess and guarantee the veracity, quality and reliability of collected data [10, 11]. There are multiple reasons for such a behaviour and here we list three common reasons:

Uncalibrated sensors: uncalibrated sensors, if ignored, can lead to severe degradation of the quality of the fields estimation [12]. Traditionally, sensors are calibrated in a controlled environment where the physical input is known and then their performance in a given calibration range is tested and verified over certain operating ranges of the environment, before such WSNs are deployed. This is infeasible for large-scale WSNs due to the prohibitive cost as well as inhomogeneity in deployment schedules. Thus, the calibration has to be done through the so-called blind or self-calibration techniques [13]. In addition, the reliability of sensor can deteriorate over time [8], making it very challenging to guarantee the quality of information even for an a priori calibrated network. Among the classical calibration models, the gain-offset response model is widely-used [14, 15]. We consider the problem of jointly estimating the calibration parameters of individual sensors as well as the spatial field values [16, 17]. 2. 2.

Compromised sensors due to malicious intent: the sensors may be physically compromised such that the sensor observations are maliciously altered to disrupt the operation of the WSN [18]. One common type of attacks is the Byzantine attack in which a hostile attacker compromises a part of the sensor network in such a manner that the Fusion Center has imperfect knowledge about whether a sensor node has been compromised [19, 20]. In such a case, there may be erroneous information incorporated into the observations from compromised sensors. 3. 3.

Unintentional misuse of sensors: in many cases data collection is done via crowd-sourcing in which private individuals install sensing stations in order to collect and share their data [21]. Since this type of data collection is not performed by professionals, in many cases, the sensors are not placed or used properly, thus introducing distortion into the measurements. Applications of crowd-sensing are becoming common recently due to the ubiquity of the Internet of Things (IoT) [22]. Those sensors could be stationary [2, 23, 24, 1, 25] or mobile [26, 27], depending on the application.

Many recent works such as [2, 28, 23, 24, 1, 25, 26] utilise the spatial-temporal correlation to reconstruct spatial physical phenomena at all locations. For example, [2] studies the placement of multi-type sensors in Gaussian spatial field to achieve optimal spatial field monitoring. To circumvent the threat to data reliability in environment monitoring systems, an estimation procedure was proposed in [28] to detect and exclude malicious sensing agents, while accurately performing spatial field reconstruction.

The main goal of this paper is to develop statistical procedures that reconstruct spatial fields using observations from sensors with possibly unknown distortions. We refer to such observations as distorted observations. We use Gaussian processes as the probabilistic model for spatial phenomena, and the sensors are assumed to follow the gain-offset distortion model with multi-modal priors to capture distortion characteristics resulted from different processes, such as natural deterioration, mis-calibration, malicious tampering, and unintentional misplacement or misuse.

The main contributions are as follows:

We develop a two-stage Bayesian inference algorithm that jointly infers the distortions of sensors and reconstructs the spatial field at all locations of interest. The algorithm estimates the distortion parameters in an empirical Bayes manner. 2. 2.

We derive the posterior distribution and the posterior predictive distribution of the model, and show that the exact computation of Bayes estimators is intractable. 3. 3.

We develop the Spatial Best Linear Unbiased Estimator (S-BLUE) for the model, which is highly computationally efficient. 4. 4.

We solve the optimization problem resulted from empirical Bayes estimation via two efficient methods, the Cross-Entropy method (CEM) and the Iterated Conditional Mode (ICM) method. 5. 5.

We analyse the computational time complexity of the proposed approaches and develop simple distributed versions of these approaches that are computationally more efficient and suitable for large-scale applications. 6. 6.

We perform synthetic data experiments as well as an experiment with real-world scenarios to validate our model and estimation procedures. The study with real-world scenarios uses a real temperature dataset from US EPA with synthetically generated distortions to show the real-world applicability of the model.

The remainder of the paper is organized as follows. We present our Bayesian sensor network model in Section II, which includes the prior distribution of the distortion parameters. In Section III, we derive the posterior distribution of the parameters as well as the posterior predictive distribution. Section IV introduces the S-BLUE and its properties. Section V introduces the approximation of the Bayes estimators via empirical Bayes, and shows that the maximization of the posterior distribution can be done through CEM and ICM. Section VI introduces the distributed approaches. In Section VII and Section VIII, we perform experiments using synthetic and real datasets. Finally, Section IX concludes the paper.

II Sensor Network Model and Assumptions

We begin by presenting the statistical model for the spatial physical phenomena, followed by the system model. The following notational convention is used throughout this paper. Boldface upper case symbols denote matrices, boldface lower case symbols denote column vectors, and standard lower case symbols denote scalars or scalar-valued functions, unless otherwise specified. All vectors are column vectors unless otherwise stated.

II-A Spatial Gaussian Random Fields Background

We model the physical phenomenon as spatially dependent continuous process with a spatial correlation structure. Such models have recently become popular due to their mathematical tractability and accuracy [1, 25, 29, 28, 3, 2]. The degree of the spatial correlation in the process increases with the decrease of the separation between two observing locations and can be accurately modeled as a Gaussian random field111We use Gaussian Process and Gaussian random field interchangeably.. A Gaussian process (GP) defines a distribution over a space of functions and it is completely specified by the equivalent of sufficient statistics for such a process, and is formally defined as follows.

Definition 1**.**

(Gaussian process [30]): Let $\mathcal{X}\subset\mathbb{R}^{d}$ be some bounded domain of a $d$-dimensional real-valued vector space. Denote by $f(\mathbf{x}):\mathcal{X}\mapsto\mathbb{R}$ a stochastic process parametrized by $\mathbf{x}\in\mathcal{X}$. Then, the random function $f(\mathbf{x})$ is a Gaussian process if all its finite dimensional distributions are Gaussian, where for any $m\in\mathbb{N}$, the random vectors $\left(f\left(\mathbf{x}_{1}\right),\ldots,f\left(\mathbf{x}_{m}\right)\right)$ has multivariate normal distribution.

We can therefore interpret a GP as formally defined by the following class of random functions:

[TABLE]

where at each point the mean of the function is $\mu(\cdot)$, and the spatial dependence between any two points is given by the covariance function (Mercer kernel) $\mathcal{C}\left(\cdot,\cdot\right)$ (see detailed discussion in [30]).

II-B Sensor Network System Model

We begin by presenting the system model followed by the prior distribution specifications.

1. A1.

Consider a random real-valued spatial phenomenon $f:\mathcal{X}\mapsto\mathbb{R}$ defined on the $d$-dimensional domain $\mathcal{X}\subset\mathbb{R}^{d}$. 2. A2.

Consider a sensor network with $N$ sensors that sense and transmit data to a Fusion Center (FC) over perfect communication channels. The spatial locations of the sensors, denoted $\left(\mathbf{x}_{n}\right)_{n=1:N}$ ($\mathbf{x}_{n}\in\mathcal{X},n=1,\ldots,N$), are known at the FC. 3. A3.

The sensor $n$ transmits $M_{n}\in\mathbb{N}$ observations to the FC. The observations $(y_{n,m})_{m=1:M_{n}}$ are generated according to the following acquisition + distortion mechanism:

[TABLE]

for $m=1,\ldots,M_{n}$, where $f\left(\mathbf{x}_{n}\right)$ is the realisation of the random field at location $\mathbf{x}_{n}$, $\epsilon_{n,m}$ represents the additive random noise at the $n$-th sensor, and $\mathcal{T}:\mathbb{R}\mapsto\mathbb{R}$ is the distortion transformation function, parametrized by ${\bm{\psi}}_{n}$. 4. A4.

The distortion transformation $\mathcal{T}$ has the following generic gain-offset form:

[TABLE]

where $a_{n}\in\mathbb{R}_{+}$ and $b_{n}\in\mathbb{R}$ represent the gain and offset of the $n$-th sensor, respectively. This gain-offset model has been widely used to describe sensor characteristics [31, 14, 15]. 5. A5.

We assume there are $K+1$ “categories” of possible distortion transformations. The auxiliary indicator random variable $Z_{n}$ indicates the category to which each sensor’s parameters ${\bm{\psi}}_{n}$ belong. We denote by $Z_{n}=0$, the default distortion transformation category, ${\bm{\psi}}_{n}={\bm{\psi}}^{0}:=(1,0)^{T}$ i.e. no distortion ($\mathcal{T}(u,{\bm{\psi}}^{0})=u$), whereas $Z_{n}=k,k\in\{1,\ldots,K\}$ indicates that the sensor $n$ belongs to the $k$-th non-default distortion transformation.

II-C Prior Distribution Specifications

• P1.

The spatial random field, $f$, is modelled as a Gaussian process (GP) $F$ with a known mean function $\mu:\mathcal{X}\mapsto\mathbb{R}$ and a known covariance function $\mathcal{C}:\mathcal{X}\times\mathcal{X}\mapsto\mathbb{R}$, that is,

[TABLE]

• P2.

The additive random noise, $\epsilon_{n,m}$, follows a normal distribution with mean zero and a fixed known variance $\varsigma^{2}$,

[TABLE]

• P3.

We place a prior distribution on $Z_{n}$, denoted $\pi({Z_{n}})$, which is a categorical distribution given by

[TABLE]

where $q^{(n)}_{k}\geq 0$ for $k=0,\ldots,K$, and $\sum_{k=0}^{K}q^{(n)}_{k}=1$.

• P4.

Each category of distortion characteristics has a distinct sub-population distribution, denoted $\pi_{k}$, which translates to the following. For $k=0,\ldots,K$,

[TABLE]

where $\pi_{0}$ is a degenerate distribution (or atom) at ${\bm{\psi}}^{0}$,

[TABLE]

For $k=1,\ldots,K$, assume that $\pi_{k}$ has density (and we slightly abuse the notation $\pi_{k}$ to also denote the density function). Thus, the prior density of ${\bm{\psi}}_{n}$ (marginalized over $Z_{n}$), denoted by $\pi({\bm{\psi}}_{n})$ is a mixture with an atom at ${\bm{\psi}}^{0}$, given by

[TABLE]

• P5.

We assume the independence among $({\bm{\psi}}_{n})_{1:N}$, and denote them collectively as ${\bm{\psi}}$. Let $\pi({{\bm{\psi}}})$ denote the prior density function of ${\bm{\psi}}$, which factorizes due to independence,

[TABLE]

The graphical structure of the proposed Bayesian model is shown in Figure 1 as a directed-acyclic-graph (DAG) using plate notations. This graphical illustration is helpful for visualizing the dependencies and conditional independence relations between parameters and random variables.

Our objective is to find an estimator $h({\bm{y}})$ for $f_{*}:=f(\mathbf{x}_{*})$, the spatial field at location $\mathbf{x}_{*}$, based on observations ${\bm{y}}:=(y_{n,m})_{n=1:N,m=1:M_{n}}$.

III Posterior of the Bayesian Model

In this section we derive the following quantities of interest, based on which the Bayesian estimators will be developed in Sections IV and V:

The posterior distribution of the model parameters ${\bm{\psi}}$, given by $p({\bm{\psi}}|{\bm{y}})$ (Theorem 1). 2. 2.

The posterior predictive distribution $p(f_{*}|{\bm{y}})$ (Theorem 2).

The following theorem gives the posterior density function of ${\bm{\psi}}$.

Theorem 1**.**

Let ${\bm{\mu}}:=\left(\mu\left(\mathbf{x}_{n}\right)\right)_{1:N}\in\mathbb{R}^{N}$ be the expected values of the $F$ process at locations $\left(\mathbf{x}_{n}\right)_{1:N}$, and let ${\bm{\mathcal{C}}}\in\mathbb{R}^{N\times N}$ be the covariance matrix, where $({\bm{\mathcal{C}}})_{ij}=\mathcal{C}\left(\mathbf{x}_{i},\mathbf{x}_{j}\right)$. For $n=1,\ldots,N$, define $g_{n}:=\sum_{m=1}^{M_{n}}y_{n,m},\;s_{n}:=\sum_{m=1}^{M_{n}}y_{n,m}^{2}.$ Let $\mathbf{M}:=\mathrm{diag}[(M_{n})_{1:N}]$, ${\bm{a}}:=(a_{n})_{1:N}$, $\mathbf{A}:=\mathrm{diag}[{\bm{a}}]$, ${\bm{b}}:=(b_{n})_{1:N}$, ${\bm{g}}:=(g_{n})_{1:N}$, ${\bm{s}}:=(s_{n})_{1:N}$. Let $\tilde{g}_{n}:=a_{n}^{-1}(M_{n}^{-1}g_{n}-b_{n})$, $\tilde{{\bm{g}}}:=(\tilde{g}_{n})_{1:N}=\mathbf{A}^{-1}(\mathbf{M}^{-1}{\bm{g}}-{\bm{b}})$, ${\bm{\Upsilon}}:={\bm{\mathcal{C}}}+\varsigma^{2}\mathbf{M}^{-1}$. Then, the log posterior density function of ${\bm{\psi}}$ is given by

[TABLE]

$p({\bm{y}})$* is a normalizing constant that is analytically intractable.*

Proof.

See Appendix A-A. ∎

Remark 1**.**

As shown in the proof of Theorem 1, the statistics $g_{n}$ and $s_{n}$ are sufficient for ${\bm{\psi}}$. In fact, as we show later, the estimators depend on ${\bm{y}}$ only through $(M_{n},g_{n},s_{n})_{n=1:N}$. Thus, from now on we take ${\bm{y}}:=(M_{n},g_{n},s_{n})_{n=1:N}$ as the summary of observations from all sensors.

The following theorem gives the posterior predictive distribution of the model.

Theorem 2**.**

The posterior predictive density is given by

[TABLE]

Let $\mu_{*}:=\mu(\mathbf{x}_{*})$, $\mathcal{C}_{*}:=\mathcal{C}(\mathbf{x}_{*},\mathbf{x}_{*})$, let ${\bm{k}}_{*}:=(\mathcal{C}\left(\mathbf{x}_{n},\mathbf{x}_{*}\right))_{n=1:N}\in\mathbb{R}^{N}$ be a column vector, and we have

[TABLE]

Proof.

See Appendix A-B. ∎

IV Spatial Best Linear Unbiased Estimator (S-BLUE)

We now derive the Spatial Best Linear Unbiased Estimator (S-BLUE). Let $l(h({\bm{y}}),f_{*})$ denote the loss function, i.e. the loss incurred when using estimator $h$ when in fact the target quantity is $f_{*}$. Let $R[\Pi,h]$ denote the Bayes risk of $h$ associated with the prior distribution $\Pi$, which is defined as the expected value of loss taken over $\Pi$, i.e.

[TABLE]

To derive the S-BLUE, we restrict the estimator to be a member of the family of linear estimators, that is $\mathcal{H}:=\{h({\bm{y}})={\bm{w}}^{T}{\bm{y}}+b\}$, where ${\bm{w}}$ is a weight vector and $b$ is an intercept, both of which do not depend on ${\bm{y}}$ or any unknown variables. Hence, the S-BLUE is defined as the optimal linear estimator under quadratic loss, $l(h({\bm{y}}),f_{*})=\left(h({\bm{y}})-f_{*}\right)^{2}$, and is given by

[TABLE]

where the expectation is taken over the joint distribution of r.v.’s $(f_{*},{\bm{y}},{\bm{\psi}})$. The next theorem shows that $\widehat{h}_{\text{S-BLUE}}$ can be expressed in closed-form.

Theorem 3**.**

$\widehat{h}_{\text{S-BLUE}}$* is given by*

[TABLE]

where $\bar{{\bm{g}}}=\mathbf{M}^{-1}{\bm{g}}$. $\mathbb{E}[\bar{{\bm{g}}}],\mathrm{Cov}[\bar{{\bm{g}}},f_{*}],\mathrm{Cov}[\bar{{\bm{g}}}]$ can all be expressed in closed-form (let $\odot$ denote matrix entry-wise multiplication),

[TABLE]

The various terms in the above equations can all be computed in closed-form, and the details are given in the proof.

Proof.

See Appendix A-C. ∎

The next corollary shows the unbiasedness property of $\widehat{h}_{\text{S-BLUE}}$.

Corollary 1**.**

$\widehat{h}_{\text{S-BLUE}}$* is unbiased, that is, $\mathbb{E}[\widehat{h}_{\text{S-BLUE}}({\bm{y}})]=\mathbb{E}[f_{*}]$.*

Proof.

It is shown via the linearity of expectation. ∎

The following corollary gives the closed-form expression of $R[\Pi,\widehat{h}_{\text{S-BLUE}}]$ (under quadratic loss).

Corollary 2**.**

Under quadratic loss, the Bayes risk associated with $\widehat{h}_{\text{S-BLUE}}$ is given by

[TABLE]

Proof.

Substituting (16) into $R[\Pi,h]=\mathbb{E}[l(h({\bm{y}}),f_{*})]$ gives the proof. ∎

The complete S-BLUE algorithm is shown in Algorithm 1. Notice that much of the computation in S-BLUE does not require ${\bm{y}}$, and thus can be performed in an “offline phase”, e.g. when the sensor network is deployed. In Algorithm 1, only Line 1 needs to be computed when the sensor measurements are taken. Below is a line-by-line analysis of the computational time complexity of Algorithm 1.

• •

Offline phase:

• –

Line 1: Evaluating $\mathbb{E}[{\bm{a}}],\mathbb{E}[{\bm{b}}]$ takes $\mathcal{O}(KN)$; evaluating $\mathbb{E}[{\bm{a}}{\bm{a}}^{T}],\mathbb{E}[{\bm{b}}{\bm{b}}^{T}],\mathbb{E}[{\bm{a}}{\bm{b}}^{T}]$ takes $\mathcal{O}(KN^{2})$.

• –

Line 1: Evaluating $\mu_{*},\mathcal{C}_{*}$ takes $\mathcal{O}(1)$; evaluating ${\bm{\mu}},{\bm{k}}_{*}$ takes $\mathcal{O}(N)$; evaluating ${\bm{\mathcal{C}}}$ takes $\mathcal{O}(N^{2})$.

• –

Line 1: Evaluating $\mathbb{E}[\bar{{\bm{g}}}],\mathrm{Cov}[\bar{{\bm{g}}},f_{*}]$ takes $\mathcal{O}(N)$; evaluating $\mathrm{Cov}[\bar{{\bm{g}}}]$ takes $\mathcal{O}(N^{2})$.

• –

Line 1: Evaluating $\widehat{{\bm{w}}}$ takes $\mathcal{O}(N^{3})$, since it involves solving a linear system; evaluating $\widehat{b}$ takes $\mathcal{O}(N)$.

• –

Line 1: Evaluating $R[\Pi,\widehat{h}_{\text{S-BLUE}}]$ takes $\mathcal{O}(N)$, since $\mathrm{Cov}[\bar{{\bm{g}}}]^{-1}\mathrm{Cov}[\bar{{\bm{g}}},f_{*}]$ has been computed in Line 1.

• •

Online phase:

• –

Line 1: Evaluating $\widehat{h}_{\text{S-BLUE}}({\bm{y}})$ takes $\mathcal{O}(N)$.

Overall, the offline phase of Algorithm 1 takes $\mathcal{O}(N^{3})$ (assuming that $K\leq N$), and its online phase takes $\mathcal{O}(N)$. It is worth noting that there are techniques to further reduce the computational complexity of the matrix inversion, e.g. through low-rank approximation (see [30]).

V Empirical Bayes Estimators

We now derive an algorithm in which we do not restrict the estimator to be linear. The idea of empirical Bayes is to plug in a point estimate $\hat{{\bm{\psi}}}$ into (12) to approximate the posterior predictive distribution, i.e. $p({\bm{\psi}}|{\bm{y}})\approx\delta_{\hat{{\bm{\psi}}}}$. This gives us the corresponding empirical Bayes estimators, which minimize the expected posterior loss, conditional on $\hat{{\bm{\psi}}}$:

[TABLE]

To complete the specification of the estimator we are required to define appropriate loss functions. We present a few widely used loss functions and their corresponding approximate Bayes estimators.

Quadratic loss function: $l_{\text{quad}}(h({\bm{y}}),f_{*})=(h({\bm{y}})-f_{*})^{2}$.

The corresponding Bayes estimator is the conditional expectation (minimum mean squared error estimator, or MMSE estimator),

[TABLE]

where $p(f_{*}|{\bm{y}})$ is given in Theorem 2. The empirical Bayes version of $\widehat{h}_{\text{MMSE}}({\bm{y}})$ is given by

[TABLE]

where $p(f_{*}|{\bm{y}},\hat{{\bm{\psi}}})$ is given in (12). 2. 2.

Absolute loss function: $l_{\text{abs}}(h({\bm{y}}),f_{*})=|h({\bm{y}})-f_{*}|$.

The corresponding Bayes estimator is the conditional median (least absolute deviation estimator, or LAD estimator),

[TABLE]

The empirical Bayes version of $\widehat{h}_{\text{LAD}}({\bm{y}})$ is given by

[TABLE] 3. 3.

$0-1$** loss function:** $l_{\text{0-1}}(h({\bm{y}}),f_{*})=\mathbbm{1}_{\{f_{*}<h({\bm{y}})\leq f_{*}+df_{*}\}}$.

The corresponding Bayes estimator is the conditional mode (maximum a posteriori estimator, or MAP estimator),

[TABLE]

The empirical Bayes version of $\widehat{h}_{\text{MAP}}({\bm{y}})$ is given by

[TABLE]

For all the aforementioned Bayes estimators, we first need to find a point estimator for ${\bm{\psi}}$. To achieve this, we find the MAP estimator of ${\bm{\psi}}$, which aims at maximizing the posterior density $p({\bm{\psi}}|{\bm{y}})$ given in Theorem 1. The MAP estimator is then given by

[TABLE]

Note that the optimization objective does not involve $p({\bm{y}})$, since it is a constant. Thus, the empirical Bayes estimators could be computed in the following two-stage algorithm:

Compute $\hat{{\bm{\psi}}}$ by solving the optimization problem $\operatorname*{arg\,max}_{{\bm{\psi}}}p\left({\bm{\psi}}|{\bm{y}}\right)$. 2. 2.

Plug in $\hat{{\bm{\psi}}}$ to compute $\widehat{h}_{\text{EB}}({\bm{y}},\hat{{\bm{\psi}}})$.

In order to solve the optimization problem in Step I, we develop two algorithms. The first approach is a stochastic optimization method named Cross-Entropy method (CEM), and the second approach is the Iterated Conditional Modes (ICM) which is based on iterative greedy search.

V-A Cross-Entropy Method (CEM)

The Cross-Entropy method (CEM) is an stochastic algorithm that is suitable for solving combinatoric or continuous optimization problems. Suppose we have a maximization problem with a unique optimizer,

[TABLE]

where $J(\cdot)$ is the objective function, $\Phi$ is the domain, and ${\bm{\varphi}}$ is the parameter vector. We solve the optimization problem by considering the level sets of the objective function $\{{\bm{\varphi}}:J({\bm{\varphi}})\geq\gamma\}$, for $\gamma\in\mathbb{R}$. When $\gamma=\hat{J}=\max_{{\bm{\varphi}}\in\Phi}J({\bm{\varphi}})$, we have $\{{\bm{\varphi}}:J({\bm{\varphi}})\geq\gamma\}=\{\hat{{\bm{\varphi}}}\}$. Next, let us define a family of probability measures $\{\mathbb{P}_{{\bm{\theta}}}:{\bm{\theta}}\in{\bm{\Theta}}\}$ on $\Phi$ with densities $\{w_{{\bm{\theta}}}:{\bm{\theta}}\in{\bm{\Theta}}\}$ that are parameterized by ${\bm{\theta}}\in{\bm{\Theta}}$. Let $\mathbb{E}_{{\bm{\theta}}}$ denote the expectation taken with respect to $\mathbb{P}_{{\bm{\theta}}}$. Let us fix ${\bm{\theta}}$ and $\gamma$, and define a rare event probability problem,

[TABLE]

Instead of approximating this probability naively by sampling from $w_{{\bm{\theta}}}$, the importance sampling method is used. Let $w_{\tilde{{\bm{\theta}}}}$ denote the importance sampler, where $\tilde{{\bm{\theta}}}\in{\bm{\Theta}}$. Importance sampling approximates the rare event probability by,

[TABLE]

where $\tilde{{\bm{\varphi}}}^{[1]},\ldots,\tilde{{\bm{\varphi}}}^{[S]}$ are $S$ independent samples generated from $w_{\tilde{{\bm{\theta}}}}$. The optimal importance sampler $w_{\hat{{\bm{\theta}}}}$ is selected through the cross-entropy criterion,

[TABLE]

where ${\bm{\varphi}}^{[1]},\ldots,{\bm{\varphi}}^{[S]}$ are $S$ independent samples generated from $w_{{\bm{\theta}}}$. Notice that the last line of (24) corresponds to the maximum likelihood estimation (MLE) of $\tilde{{\bm{\theta}}}$ when the samples are $\{{\bm{\varphi}}^{[s]}:J({\bm{\varphi}}^{[s]})\geq~{}\gamma\}$. The CEM starts from an initial sampling distribution $w_{\hat{{\bm{\theta}}}_{0}}$ and iteratively updates the threshold $\hat{\gamma}$ and the sampling distribution $w_{\hat{{\bm{\theta}}}}$. For a detailed introduction of CEM, see [32]. The complete procedure is detailed in Algorithm 2.

We can now link CEM to the MAP estimation problem in (22). We define the objective function as the (un-normalized) log-posterior conditional density,

[TABLE]

For the purpose of demonstrating the CEM, let us assume here that under prior distribution $\pi_{k}$, $(\log a_{n},b_{n})$ have a bivariate normal distribution,

[TABLE]

Note that this can be easily adapted for other prior distributions. We choose the family of sampling distributions such that,

[TABLE]

Here, we have ${\bm{\theta}}=({\bm{\theta}}^{(n)})_{n=1:N}=((r^{(n)}_{k})_{0:K}$, $(\tilde{{\bm{\nu}}}^{(n)}_{k})_{1:K}$, $(\widetilde{{\bm{\Xi}}}^{(n)}_{k})_{1:K})_{n=1:N}$. Before running the CEM algorithm, we set $\hat{{\bm{\theta}}}_{0}$ such that $w_{\hat{{\bm{\theta}}}_{0}}$ coincides with the prior distribution $\pi$. Under this setting, the optimization in Line 2 of Algorithm 2 corresponds to the MLE of ${\bm{\theta}}$, given independent samples $\{{\bm{\psi}}^{[s]}:J({\bm{\psi}}^{[s]})\geq\hat{\gamma}_{t}\}$. This decomposes into sub-problems

[TABLE]

Since $w_{{\bm{\theta}}^{(n)}}$ is a mixture distribution, the MLE does not admit a closed-form solution and we use the expectation-maximization (EM) algorithm. The EM algorithm is an iterative procedure that computes a local optimum of the likelihood function. For notational simplicity, we drop the superscripts and subscripts with $n$ for now, and denote the samples used to obtain the MLE as ${\bm{\psi}}^{[1]},\ldots,{\bm{\psi}}^{[S]}$. To apply the EM algorithm, let us first introduce the auxiliary variables. Let $z^{[s]}\in\{0,\ldots,K\}$ for $s=1,\ldots,S$ be the discrete auxiliary variables, such that

[TABLE]

This gives the marginal distributions $w_{{\bm{\theta}}}({\bm{\psi}}^{[s]})$ above. The EM algorithm starts with an initial estimate $\hat{{\bm{\theta}}}_{0}$, and iteratively updates the estimated parameter through two steps. In the expectation step (E-step), a lower bound of the log-likelihood function is constructed by first computing the conditional distributions of the auxiliary variables given the estimate of the parameters in the $t$-th iteration

[TABLE]

for $k=0,\ldots,K$, and then computing the expected value of the log-likelihood function with respect to this conditional distribution, given by

[TABLE]

In the maximization step (M-step), the estimated parameters in the $(t+1)$-th iteration are computed by maximizing the lower bound $Q({\bm{\theta}};\hat{{\bm{\theta}}}_{t})$, that is,

[TABLE]

The algorithm is summarized in Algorithm 3. For details about the EM algorithm, see [33].

To analyse the computational time complexity of CEM, let us first assume that each weighted MLE problem in Line 3 of the EM algorithm (Algorithm 3) takes $\mathcal{O}(S)$. For example, this is the case when the sampling distribution is a mixture of normal distributions. The computational time complexity of Algorithm 3 is analysed as follows:

• •

Line 3-3: Each iteration takes $\mathcal{O}(K)$, the complexity is $\mathcal{O}(SK)$.

• •

Line 3: No actual computation is performed.

• •

Line 3: Computation of $K+1$ weighted MLE takes $\mathcal{O}(SK)$.

Thus, the total complexity of the EM algorithm is $T_{\text{EM}}=\mathcal{O}(J_{\text{EM}}SK)$, where $J_{\text{EM}}$ is the number of iterations, which is usually quite small in practice. With this, we analyse the computational time complexity of the CEM estimator as follows:

• •

Preparation:

• –

Evaluation of $\log|{\bm{\Upsilon}}|$ and ${\bm{\Upsilon}}^{-1}$ takes $\mathcal{O}(N^{3})$.

• •

Algorithm 2:

• –

Line 2: Generation of $S$ samples takes $\mathcal{O}(SK)$ due to $w_{\hat{{\bm{\theta}}}_{t-1}}$ being a ($K+1$)-mixture.

• –

Line 2: Evaluation of the objective function (22) $S$ times takes $\mathcal{O}(SN^{2})$.

• –

Line 2: Computation of the sample quantile takes $\mathcal{O}(S)$.

• –

Line 2: Computation of the maximizer using Algorithm 3 takes $T_{\text{EM}}$.

• •

Inference:

• –

Evaluation of $\widehat{h}_{\text{EB}}({\bm{y}},\hat{{\bm{\psi}}})$ by (13), (14) takes $\mathcal{O}(N^{2})$.

The total complexity of the CEM estimator is $\mathcal{O}(N^{3}+J_{\text{CEM}}(SN^{2}+T_{\text{EM}}))$, where $J_{\text{CEM}}$ is the number of iterations in Algorithm 2. Since the EM algorithm converges rather quickly in practice, the complexity of CEM estimator is $\mathcal{O}(N^{3}+J_{\text{CEM}}SN^{2})$. It is worth noting that the preparation phase of the CEM estimation procedure can be run “offline”, i.e. before having access to sensor measurements. Hence, the CEM estimator has complexity $\mathcal{O}(N^{3})$ in the offline phase, and $\mathcal{O}(J_{\text{CEM}}SN^{2})$ in the online phase. The same remark on the computational complexity of matrix inversion applies here as above.

V-B Iterated Conditional Modes (ICM)

We propose a second optimization method to find the MAP estimator $\hat{{\bm{\psi}}}$ which is based on iterative greedy search. Since $\hat{{\bm{\psi}}}\in\mathbb{R}^{2N}$, the dimensionality of the optimization problem is high if $N$ is large. In addition, for $n=1,\ldots,N$, the distribution of ${\bm{\psi}}_{n}$ contains an atom. Hence, to improve the computational efficiency in these settings, we seek to reduce the complexity of the MAP estimation by reducing the global search problem to a sequence of iterative local search problems of iterated conditional modes (ICM) [34]. Let ${\bm{\psi}}_{(-n)}:=({\bm{\psi}}_{m})_{m\neq n}$. In each iteration of ICM, we fix ${\bm{\psi}}_{(-n)}$ and compute the mode of the conditional posterior distribution $\hat{{\bm{\psi}}}_{n}=\operatorname*{arg\,max}_{{\bm{\psi}}_{n}}p\left({\bm{\psi}}_{n}|{\bm{y}},{\bm{\psi}}_{(-n)}\right)$ through the conjugate gradient algorithm. ICM converges to a local maximum of the objective function in the sense that $\hat{{\bm{\psi}}_{n}}$ is the mode of the conditional posterior distribution for $n=1,\ldots,N$.

To optimize the conditional posterior distributions, we decomposed them (up to a normalizing constant) as follows, for $n=1,\ldots,N$,

[TABLE]

where ${\bm{y}}_{n}=(y_{n,m})_{m=1:M_{n}}$, ${\bm{y}}_{(-n)}=(y_{i,m})_{i\neq n,m=1:M_{i}}$. Let $\mu_{n}:=\mu\left(\mathbf{x}_{n}\right)$, ${\bm{\mu}}_{(-n)}:=(\mu_{i})_{i\neq n}$. Let $\mathcal{C}_{n}:=\mathcal{C}(\mathbf{x}_{n},\mathbf{x}_{n})$. Let ${\bm{\Upsilon}}_{(-n,n)}\in\mathbb{R}^{\left(N-1\right)}$ denote the sub-matrix of ${\bm{\Upsilon}}$ involving the cross-terms between sensor $n$ and the rest of sensors. Let ${\bm{\Upsilon}}_{(-n)}\in\mathbb{R}^{\left(N-1\right)\times\left(N-1\right)}$ denote the sub-matrix of ${\bm{\Upsilon}}$ related to sensors other than $n$. Let $\tilde{{\bm{g}}}_{(-n)}:=(\tilde{g}_{i})_{i\neq n}$. Completely analogous to Theorem 2, we have that,

[TABLE]

where

[TABLE]

One verifies that $\nu_{n}$ and $\zeta_{n}$ do not depend on ${\bm{\psi}}_{n}$. Therefore, following a derivation similar to that in Theorem 1, we have,

[TABLE]

Thus, the log-conditional likelihood as well as its partial derivatives can be efficiently evaluated. The ICM algorithm then separately treats the continuous and discrete parts of the parameter space, that is, comparing $\sup_{{\bm{\psi}}_{n}\neq{\bm{\psi}}^{0}}\log p\left({\bm{\psi}}_{n}|{\bm{y}},{\bm{\psi}}_{(-n)}\right)$ and $\log p\left({\bm{\psi}}^{0}|{\bm{y}},{\bm{\psi}}_{(-n)}\right)$.

The details of the ICM algorithm are shown in Algorithm 4. The computational time complexity of the ICM estimator is analysed as follows:

• •

Preparation:

• –

Evaluation of ${\bm{\Upsilon}}^{-1}$ takes $\mathcal{O}(N^{3})$.

• •

Algorithm 4:

• –

Line 4-4: Evaluation of $\nu_{n},\zeta_{n}$ for $n=1,\ldots,N$ takes $\mathcal{O}(N^{3})$. Notice that the computation of ${{\bm{\Upsilon}}_{(-n)}}^{-1}{{\bm{\Upsilon}}_{(-n,n)}}$ can be simplified to $\mathcal{O}(N^{2})$ via block-wise inversion once ${\bm{\Upsilon}}^{-1}$ has been computed.

• –

Line 4: Assume that the 2-dimensional optimization takes $\mathcal{O}(T_{\text{CG}})$.

• –

Line 4-4: The complexity is $\mathcal{O}(1)$.

• •

Inference:

• –

Evaluation of $\widehat{h}_{\text{EB}}({\bm{y}},\hat{{\bm{\psi}}})$ by (13), (14) takes $\mathcal{O}(N^{2})$.

The total complexity of the ICM estimator is thus $\mathcal{O}(N^{3}+J_{\text{ICM}}NT_{\text{CG}})$, where $J_{\text{ICM}}$ is the number of iterations in Algorithm 4. Similar to CEM, the preparation phase of ICM and Line 4-4 of Algorithm 4 can also be run offline. Hence, the ICM estimator has complexity $\mathcal{O}(N^{3})$ in the offline phase, and $\mathcal{O}(N^{2}+J_{\text{ICM}}NT_{\text{CG}})$ in the online phase. The same remark on the computational complexity of matrix inversion applies here as above.

To account for the multi-modality of the posterior distribution, we adopt a standard multiple start initialization strategy, that is to run ICM from a number of random initial estimates. This corresponds to running Algorithm 4 multiple times with different initial values.

VI Distributed Approaches

Now, let us consider a large scale sensor network with $I\geq 2$ clusters of sensors, where each cluster has a cluster head that locally aggregates data to be sent to the global Fusion Center. The clusters are assumed to be disjoint. We reconstruct the spatial field in a distributed manner. For a spatial location $\mathbf{x}_{*}\in\mathcal{X}$, the estimation of $f(\mathbf{x}_{*})$ is done in two steps:

Sensors within a cluster transmit the measurements to the cluster head (CH), and the CH performs a local estimation of $f(\mathbf{x}_{*})$. 2. 2.

The $I$ CHs transmit their local estimations to the FC, where local estimations are fused into the global estimation.

In this section, we develop fusion algorithms based on local S-BLUE and local empirical Bayes estimators.

VI-A Distributed S-BLUE

Suppose that each cluster head $i$ produces the local S-BLUE $\widehat{h}_{\text{S-BLUE}}^{(i)}({\bm{y}}^{(i)})$ and its Bayes risk $R[\Pi,\widehat{h}_{\text{S-BLUE}}^{(i)}]$, where ${\bm{y}}^{(i)}$ denotes the sensor measurements collected from cluster $i$. The goal is to derive a rule to fuse $\{\widehat{h}_{\text{S-BLUE}}^{(i)}({\bm{y}}^{(i)})\}_{i=1:I}$ into a single estimator $\widehat{h}_{\text{DS-BLUE}}({\bm{y}})$, where the “D” in DS-BLUE stands for “distributed”. We restrict ourselves by considering $\widehat{h}_{\text{DS-BLUE}}({\bm{y}})$ as a convex combination of $\{\widehat{h}_{\text{S-BLUE}}^{(i)}({\bm{y}}^{(i)})\}_{i=1:I}$, that is, $\widehat{h}_{\text{DS-BLUE}}({\bm{y}}):=\sum_{i=1}^{I}c_{i}\widehat{h}_{\text{S-BLUE}}^{(i)}({\bm{y}}^{(i)})$, where $c_{i}\geq 0$ for $i=1,\ldots,I$, and $\sum_{i=1}^{I}c_{i}=1$ are the constraints required to preserve the unbiasedness of DS-BLUE. The Bayes risk of $\widehat{h}_{\text{DS-BLUE}}$ is given by,

[TABLE]

where the inequality in the third line above is by the Cauchy-Schwarz inequality. The coefficients $(c_{i})_{1:I}$ are chosen to minimize the upper bound $\overline{R}[\Pi,\widehat{h}_{\text{DS-BLUE}}]$, which gives the following optimal values,

[TABLE]

It is assumed above that there is no tie among $(R[\Pi,\widehat{h}_{\text{S-BLUE}}^{(j)}])_{j=1:I}$. If there is one, the tie can be broken arbitrarily. The DS-BLUE is defined via the optimal coefficients in (28), $\widehat{h}_{\text{DS-BLUE}}({\bm{y}}):=\sum_{i=1}^{I}c^{\star}_{i}\widehat{h}_{\text{S-BLUE}}^{(i)}({\bm{y}}^{(i)})$.

The complete DS-BLUE algorithm is shown in Algorithm 5. To analyse the computational complexity of Algorithm 5, assume for now that each cluster contains at most $N_{c}$ sensor nodes. The line-by-line analysis of its computational time complexity is as follows:

• •

Offline phase:

• –

Cluster head:

Line 5: The complexity is $\mathcal{O}(N_{c}^{3})$, same as the offline phase of Algorithm 1.

• –

Fusion center:

Line 5: Evaluation of $(c_{i}^{\star})_{1:I}$ takes $\mathcal{O}(I)$.

• •

Online phase:

• –

Cluster head:

Line 5: Evaluating $\widehat{h}_{\text{S-BLUE}}^{(i)}({\bm{y}}^{(i)})$ takes $\mathcal{O}(N_{c})$.

• –

Fusion center:

Line 5: Evaluating $\widehat{h}_{\text{DS-BLUE}}({\bm{y}})$ takes $\mathcal{O}(I)$.

Overall, for cluster heads, the offline phase takes $\mathcal{O}(N_{c}^{3})$, and the online phase takes $\mathcal{O}(N_{c})$. For the fusion center, both the offline phase and the online phase take $\mathcal{O}(I)$.

VI-B Distributed Empirical Bayes Estimator

Similar to DS-BLUE, suppose that each cluster head $i$ computes an approximate posterior distribution $p(f_{*}|{\bm{y}}^{(i)},\hat{{\bm{\psi}}}^{(i)})$, where ${\bm{y}}^{(i)}$ denotes the sensor measurements collected from cluster $i$ and $\hat{{\bm{\psi}}}^{(i)}$ is a point estimate of the distortion parameters of sensors in cluster $i$. Let $\widehat{h}_{\text{EB-MMSE}}^{(i)}$ denote the empirical Bayes MMSE estimator produced by cluster head $i$. Same as DS-BLUE, let $\widehat{h}_{\text{DEB-MMSE}}({\bm{y}}):=\sum_{i=1}^{I}c_{i}\widehat{h}_{\text{EB-MMSE}}^{(i)}({\bm{y}}^{(i)})$ be a convex combination of the local estimators, where $c_{i}\geq 0$ for $i=1,\ldots,I$, and $\sum_{i=1}^{I}c_{i}=1$. We choose the fusion rule to be similar to $\widehat{h}_{\text{DS-BLUE}}$, that is,

[TABLE]

where ties are broken arbitrarily.

The complete distributed empirical Bayes algorithm is shown in Algorithm 6. For the cluster heads, the computational time complexity of this algorithm is the same as in the non-distributed version, with $N$ replaced by $N_{c}$. For the fusion center, the complexity is $\mathcal{O}(I)$.

Remark 2**.**

One downside of $\widehat{h}_{\text{DS-BLUE}}({\bm{y}})$ and $\widehat{h}_{\text{DEB-MMSE}}({\bm{y}})$ is that the reconstructed spatial field is discontinuous in space. There exists various techniques to avoid discontinuities by smoothing the reconstruction around the discontinuous boundary. However, these are left to future work.

VII Experiments with Synthetic Data

We conduct two experiments with synthetically generated data to study the performance of the methods we proposed including S-BLUE, CEM, and ICM. In Section VII-A, we study the sensitivity of the proposed methods to the strength of distortions. In Section VII-B, we perform a realistic simulation and analyse the overall performance of the proposed methods. In the studies, the two empirical Bayes-based estimators (CEM and ICM) use the quadratic loss function and hence correspond to $\widehat{h}_{\text{EB-MMSE}}$. The proposed methods are compared to two baselines, the “oracle” case in which the distortions ${\bm{\psi}}$ are known exactly, and the ”naive” case in which distortions are disregarded in the prediction, i.e. $\hat{{\bm{\psi}}}_{n}={\bm{\psi}}^{0}$, for $n=1,\ldots,N$.

VII-A Synthetic Experiment 1: Homogeneous Distortion Characteristics

In this experiment, we study an ideal scenario where the signal-to-noise ratio (SNR) is high and observations are plentiful. We fix the number of observations per sensor to be 50, and the SNR to be 15dB. Notice that for the ease of comparison, all SNRs are measured at the sensor level, that is, the SNRs of aggregated observations. Under the i.i.d. noise assumption, we define $\text{SNR}=10\log_{10}\left(\frac{M\mathrm{Var}(F)}{\varsigma^{2}}\right)$, where $\mathrm{Var}(F)$ denotes the signal variance, $\varsigma^{2}$ denotes observation noise variance, and $M$ denotes the number of observations per sensor.

We simulate a spatial field defined on the two-dimensional square $\mathcal{X}=[0,1]^{2}$ with mean 10, i.e. $\forall\mathbf{x}\in\mathcal{X}$, $\mu(\mathbf{x})=10$, and a Matérn covariance function with $\nu=3/2$, $\mathrm{Var}(F)=100$, and length scale=0.3, i.e. $\forall\mathbf{x},\mathbf{x}^{\prime}\in\mathcal{X}$, $\mathcal{C}(\mathbf{x},\mathbf{x}^{\prime})=100\left(1+\frac{\sqrt{3}\|\mathbf{x}-\mathbf{x}^{\prime}\|}{0.3}\right)\exp\left(-\frac{\sqrt{3}\|\mathbf{x}-\mathbf{x}^{\prime}\|}{0.3}\right).$ Here $\|\mathbf{x}-\mathbf{x}^{\prime}\|$ corresponds to the Euclidean distance. The contour plot of the simulated spatial field is shown in Figure 2(a).

Subsequently, 100 sensors are randomly placed in the square. 50 out of the 100 sensors are fixed to have identical distortion parameters, and the rest are set to have the default transformation parameters ${\bm{\psi}}^{0}$, that is, non-distorting. For the sensors with distortions, we first fix the offset parameter $b_{n}$ at 5, and vary the gain parameter $a_{n}$ from 1 to 1.6. Then we fix the gain parameter at 1.2, and vary the offset parameter from 0 to 12. With each setting of distortion parameters, 100 sets of noisy observations are randomly simulated. For each set of observations, the three proposed methods: S-BLUE, CEM, and ICM, along with the two baselines oracle and naive, are used to reconstruct the spatial field. Here, weakly informative prior for the distortion parameters is used, which has a single category ($K=1$) given by $q^{(n)}_{1}=0.5$ for all $n$, where under $\pi_{1}$, $a_{n}\sim\log\mathcal{N}(0.25,0.1^{2})$, $b_{n}\sim\mathcal{N}(6,3^{2})$. The reconstruction accuracy is evaluated by the mean-squared-errors (MSE) at a $100\times 100$ grid on $[0,1]^{2}$.

Figure 3 shows the reconstruction accuracy averaged over 100 realizations. For better interpretability, the ratio between the MSE and the prior variance, referred to as the relative MSE, is shown. Error bars in Figure 3 indicate the 95% Student’s $t$-confidence interval of the relative MSE estimated from the 100 realizations. Error bars in all subsequent figures indicate the 95% Student’s $t$-confidence interval of the respective underlying quantity. From Figure 3, one observes that the MSE of the oracle stayed constant, as expected, while the MSE of the naive baseline increased rapidly when the distortion parameters increased. The MSE of S-BLUE first decreased and then increased slightly. This is due to the way the prior distribution of distortions were set up. The prior mean of $a_{n}$ and $b_{n}$ were 1.29 and 6.0, respectively. Since S-BLUE makes predictions based purely on the prior information, its performance is best when the actual distortion is closest to the prior mean. The two empirical Bayes-based methods showed decreasing MSE when the gain parameter increased and slightly increasing MSE when the offset parameter increased. The reason is that since the gain parameter affects both the location and the spread of the observations, while the offset parameter only affects the location, the gain in the distortion was more noticeable and thus easier to detect. The MSEs of CEM and ICM were almost identical. Figure 4 shows the average false positive rate (FPR) and false negative rate (FNR) of CEM and ICM. The FPR is defined as the proportion of non-distorting sensors that were estimated to be distorting, and the FNR is defined as the proportion of distorting sensors that were estimated to be non-distorting. The FPR and FNR of the two methods were almost identical. Notice in addition that the FNR was high when the gain was 1 and the offset was 5. This was caused by the short length scale (0.3), which made detection of the offset hard due to the low spatial correlation. Finally, all error bars are narrow, indicating that difference between the performance of different methods are statistically significant. To further confirm this, we show the maximum absolute deviation from the relative MSE of the five methods in the second column of Table I. One checks that the deviations are small, indicating that the performance is stable across realizations.

VII-B Synthetic Experiment 2: Inhomogeneous Distortion Characteristics

In the second synthetic experiment, we study a realistic scenario where each sensor has different distortion parameters, and vary the SNR as well as the number of observations.

We again simulate a spatial field defined on the two-dimensional square $[0,1]^{2}$. This time, however, the length scale is set to be 0.5, and the spatial correlation decays at a slower rate. The contour plot of the simulated spatial field is shown in Figure 2(b). The 100 sensors are placed at the same locations as in the synthetic experiment 1. 50 out of the 100 sensors are randomly selected to have the different distortion parameters generated from the following prior with three categories ($K=3$), given by $q^{(n)}_{1}=q^{(n)}_{2}=q^{(n)}_{3}=\frac{1}{6}$ for all $n$, where under $\pi_{1}$, $a_{n}\sim\log\mathcal{N}(-0.4,0.05^{2})$, $b_{n}\sim\mathcal{N}(0,0.2^{2})$, under $\pi_{2}$, $a_{n}\sim\log\mathcal{N}(0.2,0.05^{2})$, $b_{n}\sim\mathcal{N}(0,0.2^{2})$, under $\pi_{3}$, $a_{n}\sim\log\mathcal{N}(0,0.05^{2})$, $b_{n}\sim\mathcal{N}(10,2^{2})$, and the rest of the sensors are set to be non-distorting. The distortion parameters are generated and fixed in this experiment. After that, we randomly simulate 100 sets of noisy observations. For each set of observations, we test the proposed methods along with the baselines as in the synthetic experiment 1.

Figure 5 shows the relative MSE averaged over 100 realizations, with different number of observations per sensor and different SNR. Figure 6 shows the FPR and FNR of CEM and ICM averaged over 100 realizations. Observe that the MSE of the baselines and S-BLUE did not change with different number of observations per sensor because they depend only on the mean of observations from each sensor. CEM and ICM, on the other hand, benefited from having access to more observations. The naive baseline showed a peculiar trend that first decreased and then increased when SNR increased. The reason is that the naive estimator is a linear combination of the prior mean and the observations, where the weights depend on the SNR. With high SNR, the naive estimator placed a high weight on the distorted observations, thus making the MSE high. The MSE of S-BLUE decreased steadily when SNR increased, and eventually flattened out. The MSE of CEM and ICM depended highly on the number of observations. CEM had the best reconstruction quality among the proposed methods when observations were plentiful or when the SNR was high. In comparison, ICM performed considerably worse than CEM when the number of observations was 5 and 20. This was due to the higher spatial correlation, which made the dependency between distortion parameters higher in the posterior and reduced the effectiveness of iterative greedy search. Notice that the FPR of CEM and ICM was close to 1 when the SNR was low and the observations were scarce. This indicates that CEM and ICM estimated all of the sensors as distorting. The reason was that with high noise variance, the likelihood had a flat shape, and thus the posterior mode did not contain non-distorting sensors. It can also be observed that the FPR and FNR were low when the number of observations was 100 and the SNR was high. Same as in the synthetic experiment 1, the error bars are narrow, indicating the statistical significance of the differences. The third column of Table I shows the maximum absolute deviation from the relative MSE in the synthetic experiment 2. The deviations are larger compared to the synthetic experiment 1 due to the case with small number of samples and low SNR. Nonetheless, this shows that the performance of the methods is stable across realizations. To further demonstrate the performance of the proposed methods under different distortion parameters, we repeat the above experiment. This time, instead of fixing the distortion parameters across realizations, the distortion parameters are independently randomly generated in each of the 100 realizations. The result is shown in Figure 7. Since Figure 5 and Figure 7 look very similar, we confirm that the performance is stable across a range of distortion parameters.

Finally, to show the effect of the proportion of distorting sensors, an additional experiment is performed. Figure 8 shows the relative MSE of S-BLUE, CEM, ICM and the baselines when the number of observations is fixed at 100, the SNR is fixed at 20dB, and the proportion of distorting sensors varies from 0 to 1. These results were averaged across 100 independent realizations of distortion parameters and random noises. From Figure 8, it can be observed that when no sensor was distorting, the relative MSE of all methods coincided. However, as the proportion of distorting sensors increased, the relative MSE of the naive baseline increased rapidly, the relative MSE of S-BLUE increased slowly, the relative MSE of CEM and ICM increased only slightly and flattened out eventually, and the oracle, as expected, was unaffected by the change of proportion. The reason is that increasing the proportion of distorting sensors resulted in an increase in prior model uncertainty that negatively affected S-BLUE, while CEM and ICM used the information from the observations to estimate the distortion parameters and thus were only slightly affected.

VIII Experiments with Real Data

We study the 2017 US temperature dataset from the US EPA222https://aqs.epa.gov/aqsweb/airdata/download_files.html, retrieved on 5 June 2018.. The dataset contains 309,226 rows, with the following fields:

• •

State.Code: the numerical code of the state.

• •

County.Code: the numerical code of the county.

• •

Site.Num: the numerical code of the monitoring site.

• •

Longitude: the longitude of the site.

• •

Latitude: the latitude of the site.

• •

Date.Local: the local date on which the temperature measurement was taken.

• •

X1st.Max.Value: the maximum hourly temperature measurement of a day.

VIII-A Preprocessing

The first step of preprocessing is to remove the irrelevant fields from the dataset. For temperature measurements, we take the maximum hourly measurement of every day. The measurements are converted from Fahrenheit into Celsius. We noticed that there is an obvious outlier which corresponds to 125$\degree$C (the next highest measurement is 55$\degree$C) with state code 6, county code 79, site number 5 on day 161. We also noticed another potential outlier corresponding to $-17.89\degree$C ($0\degree$F) in June, while the next lowest temperature from May to July is 1.1$\degree$C. This measurement has state code 38, county code 93, site number 101 on day 181. Therefore, these two measurements are removed from the dataset.

In the dataset, there are 830 monitoring sites in total. Figure 9 is a scatter plot showing the spatial locations of these monitoring sites. Not all monitoring sites have taken measurements on all 365 days. We refer to these as missing measurements. Out of the 830 monitoring sites, 55 contain more than 30% of missing measurements. Overall, 16,980 measurements are missing, which is 5.6% of the 302,950 measurement (365$\times$830).

Remark 3**.**

We have noticed that some monitoring sites have reported temperature measurements by different types of instrument, which explains why we have slightly fewer measurements after we sort measurements into the site & date format. However, measurements on the same day at the same monitoring site with different instruments tend to be close hence we choose an arbitrary set whenever a monitoring site reports multiple sets of measurements.

Figure 10 shows the histogram and the normal quantile-quantile (Q-Q) plot of the measurements. From Figure 10(a), one can see that the measurements contain quantization artifacts. Nonetheless, the measurements are treated as continuous. From Figure 10(b), one sees that the dataset does not contain abnormally large or small values, and is slightly left-skewed compared to a normal distribution. Overall, the normality assumption holds approximately, as we have not yet considered the seasonal shift of temperature. Table II shows the summary statistics of these measurements.

VIII-B Smoothing of Spatial Field

The first step of the spatial analysis is to model the daily observations as noisy samples from a Gaussian process without distortions and estimate its hyperparameters. For subsequent analyses, we take only data from the 20 days between day 181 to day 200 from the dataset, and restrict ourselves to the contiguous United States (that is, excluding Hawaii and Alaska) since Hawaii and Alaska are far away from the rest of the Unites States. As for the covariance function, we again choose to use the Matérn covariance function with $\nu=3/2$ to allow for high flexibility while keeping the spatial field mean-square differentiable. For each calendar day, we estimate the signal mean, signal variance, noise variance and length-scale via maximum marginal likelihood estimation (see Chapter 5 of [30] for details) using all available observations on that day. Then, the median of the estimated values on 20 days are taken as the estimated hyperparameters of the GP. The estimated hyperparameters are shown in Table III.

Using the estimated hyperparameters, we reconstruct the spatial field for each calendar day at both a $100\times 100$ grid of locations and the 824 sensor locations by computing the posterior mean and the posterior covariance matrix, conditional on all the available observations on that day. Since the posterior mean is usually much smoother compared to the actual spatial field, a sample spatial field is generated from the posterior distribution to make it realistic. The generated spatial field at a grid of locations is treated as the ground truth of the spatial field, and the generated field intensities at the 824 sensor locations are treated as as the noise-free sensor reading from which noisy and distorted observations are generated.

VIII-C Experimental Settings

In this experiment, in addition to the five methods tested in the synthetic experiments, we evaluate the performance of the distributed version of S-BLUE, CEM, and ICM introduced in Section VI. To do so, we first divide the 824 sensor locations into 8 disjoint clusters via hierarchical clustering with great-circle distance and complete linkage (see e.g. Section 14.3.12 of [35]). The 8 resulting clusters are indicated in Figure 9 using 8 different colors. Subsequently, we apply the distributed approaches to this partitioned sensor network.

In order to evaluate the reconstruction accuracy of the proposed methods, $25\%$ of the sensor locations are left out as the test set. For the 618 remaining locations, we randomly generate the distortion parameters from three different settings. In the three settings, each sensor has a respective probability of 0.3, 0.5 and 0.7 to introduce distortion. Under all three settings, the distorting sensors have one of the three following categories, with equal probabilities. Under $\pi_{1}$, $a\sim\log\mathcal{N}(-0.4,0.05^{2})$, $b\sim\mathcal{N}(0,0.2^{2})$, under $\pi_{2}$, $a\sim\log\mathcal{N}(0.2,0.05^{2})$, $b\sim\mathcal{N}(0,0.2^{2})$, under $\pi_{3}$, $a\sim\log\mathcal{N}(0,0.05^{2})$, $b\sim\mathcal{N}(10,2^{2})$. For each of the three settings, we first randomly generate the distortion parameters for each sensor. Then we randomly simulate noisy observations at each sensor location and subsequently apply the corresponding distortions. In addition, we examine the effect of the number of observations and effective SNR. Specifically, we examine the four following cases:

10 observations per sensor, SNR=5dB; 2. 2.

10 observations per sensor, SNR=15dB; 3. 3.

50 observations per sensor, SNR=5dB; 4. 4.

50 observations per sensor, SNR=15dB.

VIII-D Results and Discussion

Figure 11 shows the relative MSE averaged over 20 days, and Figure 12 shows the FPR and FNR of CEM and ICM, under each setting. Similar to the synthetic experiment 2, we observe that the MSE of the baselines and S-BLUE did not change with different number of observations per sensor. CEM and ICM, on the other hand, benefited from having access to more observations.

The naive baseline showed worse MSE when SNR was high. The reason was the same as in the synthetic experiment 2. In addition, the MSE of naive baseline is clearly affected by a higher proportion of distorting sensors.

S-BLUE showed stable accuracies across all settings. Since S-BLUE does not estimate the distortion parameters, the error mainly resulted from smoothing. This can be seen clearly from the reconstructed spatial fields in Figure 13, which will be discussed later.

Compared to the simple method S-BLUE, the more sophisticated methods CEM and ICM have the additional benefit of being able to estimate the distortion parameters. CEM consistently outperformed ICM, which is consistent with what we observed in Section VII-B. This possibly indicates the ineffectiveness of the iterative greedy search method with high-dimensional mixed discrete-continuous optimization problems. While CEM benefited slightly from higher SNR, the performance of ICM deteriorated with higher SNR when there were 10 observations per sensor, though they both benefited from more plentiful observations. Looking at Figure 12, it can be observed that ICM had significantly higher FPR and slightly lower FNR compared to CEM. More observations and higher SNR have the effect of decreasing the FPR of both CEM and ICM.

Remark 4**.**

We noticed that if the SNR is set to be higher than 20dB, the performance of both CEM and ICM deteriorates greatly. With a high SNR (hence low noise variance), the posterior density of the model tends to be highly irregular. Thus the optimization problem might be highly ill-posed. In general, there is no universal solution to this. Some form of relaxation might help to improve the performance, but that is not examined in this work.

In the low SNR settings, the distributed approaches performed slightly worse compared to their centralized counterparts. In the high SNR settings, the distributed S-BLUE performed slightly worse compared to the centralized S-BLUE. However, the distributed CEM and ICM performed better compared to their centralized version. This indicates that both CEM and ICM suffered from convergence issues when dealing with the high-dimensional optimization problem. In the distributed version, the problem is decomposed into sub-problems of lower dimensionality, which are easier to solve. This experiment shows that the distributed approaches we proposed are the ideal candidates for solving such spatial field reconstruction problem in large-scale settings.

Table IV shows the maximum absolute deviation from the relative MSE in this experiment. It shows that S-BLUE and the distributed CEM had stable performance throughout the 20 days. ICM and the centralized CEM are less stable, presumably due to the difficulty in the optimization procedures.

Finally, let us examine the reconstructed spatial fields. Figure 13 shows the heat maps of spatial fields reconstructed by centralized and distributed versions of S-BLUE, CEM, ICM, and the two baselines with the settings: proportion of distorting sensors is 0.7, 10 observations per sensor, SNR=15dB. The ground truth spatial field is also included for reference. First, notice that the ground truth contained much more details compared to the reconstructions, because all of the estimators here have the smoothing effect. Overall, the reconstruction of S-BLUE is much smoother and contains fewer details. This is due to nature of S-BLUE, as it does not estimate the distortion parameters. CEM and ICM, on the other hand, preserved many of the details, and their reconstructions were overall close to the one produced by the oracle. The naive method, however, produced a noticeably inaccurate reconstruction. As previously mentioned, the reconstructions of the distributed approaches are discontinuous and the discontinuity can be observed at the boundary of the clusters. Nonetheless, they produced accurate reconstructions of the spatial field.

Notice that the above analysis used only 20 consecutive days of data. When a longer time period is considered, the seasonal variation of the underlying spatial field must be taken into consideration. To demonstrate this, we use the data from each of the twelve months to estimate the GP hyperparameters by the same approach as above. Subsequently, the spatial field on the last day of each month is reconstructed by centralized CEM, using noisy and distorted sensor readings at the 618 locations generated in the same way as in the above experiment. These reconstructions along with the ground truths are shown in Figure 14. It can be seen that there is an obvious seasonal effect on the temperature throughout the year. Notice that the reconstructions closely resemble the ground truths, which indicates that CEM worked well when the characteristics of the underlying spatial field varied throughout the year.

IX Conclusion

This paper addressed the problem of spatial field reconstruction based on distorted sensor readings. A new spatial field model based on a mixture of Gaussian process experts was developed. We developed two approaches to solve the inference problem. The first approach uses a linear Bayes estimator named the Spatial Best Linear Unbiased Estimator (S-BLUE), which is a low-complexity algorithm relying only on prior information. The second approach is a two-stage algorithm based on empirical Bayes, in which the unknown distortion parameters of the sensors are estimated based on distorted observations. We developed two optimization procedures for the two-stage algorithm, the first one is based on the Cross-Entropy method (CEM) and the second one is based on the Iterated Conditional Mode (ICM) which is an iterative greedy search procedure. In addition, the distributed versions of S-BLUE and empirical Bayes estimators were developed to improve the computational efficiency in large-scale applications. We preformed two synthetic experiments as well as an experiment based on real temperature data from US EPA with synthetically generated distortions to assess the spatial field reconstruction accuracy of the proposed approaches. The results showed significant improvement compared to the estimation approach that neglects sensor distortions.

Acknowledgments

The research was conducted under the Cooling Singapore project, funded by Singapore’s National Research Foundation (NRF) under its Virtual Singapore programme.

Appendix A Proof of Theorems

A-A Proof of Theorem 1

Proof.

For $n=1,\ldots,N$, the log model likelihood given ${\bm{f}}:=(f(\mathbf{x}_{n}))_{n=1:N}$ is given by

[TABLE]

which depends on $y_{n,1:M_{n}}$ only through the statistics $g_{n}$ and $s_{n}$. Therefore, $({\bm{g}},{\bm{s}})$ are sufficient for $({\bm{f}},{\bm{\psi}})$. The joint density of $({\bm{y}},{\bm{f}})$ conditional on ${\bm{\psi}}$ is given by

[TABLE]

where ${\bm{\gamma}}=\varsigma^{-2}\mathbf{M}\tilde{{\bm{g}}}+{\bm{\mathcal{C}}}^{-1}{\bm{\mu}},\mathbf{Z}=\varsigma^{-2}\mathbf{M}+{\bm{\mathcal{C}}}^{-1}$. Integrating over ${\bm{f}}$, we deduce that,

[TABLE]

We have by Woodbury’s formula that,

[TABLE]

We also have by ${\bm{\gamma}}=\varsigma^{-2}\mathbf{M}\tilde{{\bm{g}}}+{\bm{\mathcal{C}}}^{-1}{\bm{\mu}}$ that,

[TABLE]

Substituting (32) into the first term on the right-hand side of (35), using the fact that $\mathbf{M}\mathbf{Z}^{-1}{\bm{\mathcal{C}}}^{-1}=({\bm{\mathcal{C}}}\mathbf{Z}\mathbf{M}^{-1})^{-1}=(\varsigma^{-2}{\bm{\mathcal{C}}}+\mathbf{M}^{-1})^{-1}=\varsigma^{2}{\bm{\Upsilon}}^{-1}$ for the second term of (35), substituting (33) into the third term of (35), and finally substituting the resulting formula as well as (34) into (31), one gets the following after simplification (recall that $\tilde{{\bm{g}}}:=\mathbf{A}^{-1}(\mathbf{M}^{-1}{\bm{g}}-{\bm{b}})$, ${\bm{\Upsilon}}:={\bm{\mathcal{C}}}+\varsigma^{2}\mathbf{M}^{-1}$),

[TABLE]

Notice that after integrating out ${\bm{f}}$, $({\bm{g}},{\bm{s}})$ are still sufficient for ${\bm{\psi}}$.

By Bayes’ rule, $p({\bm{\psi}}|{\bm{y}})=\frac{p({\bm{y}}|{\bm{\psi}})\pi({\bm{\psi}})}{p({\bm{y}})}$, where $p({\bm{y}})=\int p({\bm{y}}|{\bm{\psi}})\pi({\bm{\psi}})d{\bm{\psi}}$ is the normalizing constant that is analytically intractable. The proof is now complete. ∎

A-B Proof of Theorem 2

From (30) we see that $\log p({\bm{y}},{\bm{f}}|{\bm{\psi}})$ has the following form,

[TABLE]

where $q_{1},q_{2},q_{3}$ are some functions of the corresponding parameters. Similarly, we can deduce that $\log p({\bm{y}},{\bm{f}},f_{*}|{\bm{\psi}})$ has the following form,

[TABLE]

Therefore,

[TABLE]

Thus, we have deduced that $p(f_{*}|{\bm{y}},{\bm{\psi}})$ depends on ${\bm{y}}$ only through ${\bm{g}}$ (or equivalently, $\tilde{{\bm{g}}}$), i.e. $p(f_{*}|{\bm{y}},{\bm{\psi}})=p(f_{*}|{\bm{g}},{\bm{\psi}})=p(f_{*}|\tilde{{\bm{g}}},{\bm{\psi}})$. We also have that conditional on ${\bm{\psi}}$, $(\tilde{{\bm{g}}},f_{*})$ are jointly Gaussian, with $\mathbb{E}[f_{*}|{\bm{\psi}}]=\mu_{*}$, $\mathbb{E}[\tilde{{\bm{g}}}|{\bm{\psi}}]={\bm{\mu}}$, $\mathrm{Cov}[f_{*}|{\bm{\psi}}]=\mathcal{C}_{*}$, $\mathrm{Cov}[\tilde{{\bm{g}}},f_{*}|{\bm{\psi}}]={\bm{k}}_{*}$, $\mathrm{Cov}[\tilde{{\bm{g}}}|{\bm{\psi}}]={\bm{\Upsilon}}$, which can be easily verified. Thus, we have that $(f_{*}|{\bm{y}},{\bm{\psi}})\sim\mathcal{N}(\bar{f}_{*},\sigma^{2}_{*})$, where $\bar{f}_{*}=\mu_{*}+{\bm{k}}_{*}^{T}{\bm{\Upsilon}}^{-1}\left(\tilde{{\bm{g}}}-{\bm{\mu}}\right)$, $\sigma^{2}_{*}=\mathcal{C}_{*}-{\bm{k}}_{*}^{T}{\bm{\Upsilon}}^{-1}{\bm{k}}_{*}$. The distribution of $(f_{*},{\bm{\psi}}|{\bm{y}})$ is given by $p(f_{*},{\bm{\psi}}|{\bm{y}})=p(f_{*}|{\bm{y}},{\bm{\psi}})p({\bm{\psi}}|{\bm{y}})$, and the posterior predictive distribution $(f_{*}|{\bm{y}})$ is obtained by marginalizing over ${\bm{\psi}}$, $p(f_{*}|{\bm{y}})=\int p(f_{*}|{\bm{y}},{\bm{\psi}})p({\bm{\psi}}|{\bm{y}})d{\bm{\psi}}$. The proof is now complete.

A-C Proof of Theorem 3

Proof.

First, we claim that any linear estimator must have the form $h({\bm{y}})={\bm{w}}^{T}\bar{{\bm{g}}}+b$, where ${\bm{w}}\in\mathbb{R}^{N},b\in\mathbb{R}$. Due to the linearity of $h$ in observations $(y_{n,m})_{n=1:N,m=1:M_{n}}$, ${\bm{s}}$ is not involved. For $n=1,\ldots,N$, the weights of $(y_{n,m})_{m=1:M_{n}}$ must be the same due to symmetry. This proves the claim.

Under quadratic loss, for $h\in\mathcal{H}$, $R[\Pi,h]$ is given by

[TABLE]

Differentiating $R[\Pi,h]$ with respect to ${\bm{w}}$ and $b$, we get

[TABLE]

Setting the partial derivatives to $\mathbf{0}$ and solving the equations gives the optimal weight vector and intercept,

[TABLE]

One can verify that $(\widehat{{\bm{w}}},\widehat{b})$ is indeed minimizing the Bayes risk. Hence,

[TABLE]

The terms $\mathbb{E}[f_{*}],\mathbb{E}[\bar{{\bm{g}}}],\mathrm{Cov}[\bar{{\bm{g}}},f_{*}],\mathrm{Cov}[\bar{{\bm{g}}}]$ can all be expressed in closed-form. The closed-form expressions and the details of the computation are given below. Let $\odot$ denote matrix entry-wise product.

[TABLE]

In the terms $\mathbb{E}[{\bm{a}}]$, $\mathbb{E}[{\bm{b}}]$, $\mathbb{E}[{\bm{a}}{\bm{a}}^{T}]$, $\mathbb{E}[{\bm{b}}{\bm{b}}^{T}]$, $\mathbb{E}[{\bm{a}}{\bm{b}}^{T}]$, the expectations are evaluated entry-wise. For example, entries of $\mathbb{E}[{\bm{a}}]$ are given by

[TABLE]

Entries of $\mathbb{E}[{\bm{a}}{\bm{a}}^{T}]$ are given by

[TABLE]

Entries of $\mathbb{E}[{\bm{b}}]$, $\mathbb{E}[{\bm{b}}{\bm{b}}^{T}]$ and $\mathbb{E}[{\bm{a}}{\bm{b}}^{T}]$ can be evaluated similarly. With the above equations, we are able to evaluate $\widehat{h}_{\text{S-BLUE}}$ efficiently. ∎