Strong consistency of a kernel-based rule for spatially dependent data Open Access

In many applications one needs to classify spatial data that have been collected incompletely. The classification of incomplete-data problem, in which certain features are missing from particular feature vectors, exists in a wide range of fields, including image labeling, computer vision and others. For example, in the remote sensing technology, because of the internal malfunction of satellite sensors and poor atmospheric conditions such as thick cloud, the acquired remote sensing images often suffer from missing information at certain pixels and one wants to classify these pixels using the information in the nearest identified pixels. Many existing classification algorithms assume either certain parametric distributions for the data or certain forms of separating curves or surfaces. These parametric classifiers are suboptimal and of limited use in practical applications where little information about the underlying distributions is available a priori. In comparison, nonparametric classifiers are usually more flexible in accommodating different data structures, and are hence more desirable. [21] has proposed a nonparametric approach allowing to include contextual features for classifying missing spatial data and has investigated the consistency of the classifier under mild conditions. In nonparametric spatial estimation, the existing works concern mainly the estimation of a probability density and regression functions, see the key references: [2–4,15] and [14]. More recently, [5] has proposed a kernel spatial density estimator allowing for the analysis of spatial clustering. In this work, we establish strong consistency of the classifier proposed by [21] and then, we check its performance with simulation studies and applications. We consider a strictly stationary random field ${(Xi,Yi)}i∈ℤN$ defined on some probability space $(Ω,ℱ,ℙ)$ and taking values in $ℝd×{0,…,M}$⁠, for some integer $M≥1$⁠. In the problem of classification, for each $i∈ℤN$⁠, $Xi$ is a vector of features and $Yi$ is the label (class) of $Xi$⁠. A point $i=(i1,…,iN)∈ℤN$ will be referred to as a site. For $n=(n1,…,nN)∈(ℕ∗)N$⁠, we define the rectangular region $In$ by $In={i∈ℤN:1≤ik≤nk,∀k=1,…,N}$⁠. We will write $n→∞$ if $mink=1,…,Nnk→∞$⁠. Define $n^=n1×⋯×nN=card(In)$ and assume that the random field is observed on a subset $Sn⊂In$ with $In−Sn$ is a bounded set for $n^$ large enough. When processing a particular site, its features are not used at all, but only the features of its neighbors will be considered. In other words, we wish to predict the label $Yj$ of a new site $j$ based only on observations in a vicinity, say $νj⊂Sn$⁠, where the set $νj$ is not containing $j$⁠. Let $νj=j+ν$⁠, where $ν⊂ℤN$ is a fixed bounded set of sites not containing $0$ with $card(ν)=l$ (⁠$l$ is also the cardinal of each $νj$⁠). We assume that $X(j)={Xi:i∈νj}$ is a random vector taking values in $ℝd˜$ with $d˜=ld$⁠, and that the components of $X(j)$ are ordered according to an arbitrary order on indices, for example the lexicographic order. The pair $(X(j),Yj)$ may be completely described by $μ$⁠, the probability measure for $X(j)$⁠, and $η(x)$⁠, the regression of $Yj$ on $X(j)=x$⁠. Assume that for each $i∈ℤN$⁠, $(X(i),Yi)$ has the same distribution as the pair $(X(1),Y1)$⁠. We will create a classifier $g:ℝd˜→{0,…,M}$ mapping $X(j)$ into the predicted label of $Xj$⁠. The error rate, or risk, of a rule $g$ is $L(g)=ℙ{g(X(j))≠Yj}$⁠. This is minimized by the rule

whose error rate $L∗=L(g∗)$ is called the Bayes-optimal risk and $g∗(x)$ is called the Bayes rule. Clearly, $g∗(x)$ predicts the label $Yj$ of the site $j$ using only $x$⁠, the value of $X(j)$⁠, while the features vector $Xj$ does not affect the classification procedure at all. This means that $g∗(x)$ well work event if $Xj$ is completely missing. Unfortunately, we cannot use (1.1) directly because it depends on the distribution of $(X(j),Yj)$ which is generally unknown. So, we take $Jn={i∈Sn:νi⊂Sn}$ and we use the training data $Dn={(Xi,Yi):i∈Jn}$ to construct a classifier $gn(x)$⁠. We consider the classifier $gn(x)$ obtained by extending the classifier of [21] to the multi-class case as follows:

where $1A$ denotes the indicator of the set $A$⁠, the kernel $K:ℝd˜→ℝ+$ is a density function on $ℝd˜$⁠, and $bn$ is a sequence of bandwidths tending to zero as $n$ tends to infinity. In one hand, the sum in (1.2) is taken over $Jn$ instead of $Sn$ just to ensure that $X(i)$ always exists and that the sums make sense. On the other hand, for each new site $j∉Sn$⁠, the classifier $gn(x)$ predicts the missing label $Yj$ independently of its features vector $Xj$ which does not belong neither to the training sample $Dn$ nor to the components set of $X(j)$⁠. Consequently, $gn(x)$ may classify $j$ even if its own features vector $Xj$ is completely missing and that makes our method exhibit good performance in comparison with the classical spatial Markovian model. [6] proposes a nonparametric approach to extend the result of [2] to the non-Markovian case by using two kernels in the estimator in order to control both the distance between observations and that between spatial locations without using a specific vicinity for the non-observed site. This latter approach may be developed to classify spatial data but it does not work when one wants to classify sites with missing or incomplete features. Let $Ln=ℙ{gn(X(j))≠Yj|Dn}$ be the error probability of $gn(x)$⁠. Generally, we cannot hope to design a classifier that achieve the Bayes error probability $L∗$ but it is possible that the limit behavior of $Ln$ compares favorably to $L∗$⁠. This idea is encapsulated in the notion of consistency.

Let $(Ω,ℱ,ℙ)$ be a probability space and let $A$ and $B$ be two sub $σ$-fields of $ℱ$⁠. The $α$-mixing coefficient between $A$ and $B$ is defined by

and the $β$-mixing coefficient is defined by

Let $(Zi)i∈ℤN$ be a random field on (Ω$,ℱ,ℙ)$ and taking values in some space (Ω^′$,ℱ′$⁠).

This section is a collection of technical lemmas which will be used to prove the strong consistency result stated in Theorem 4.1. Let $‖.‖r$ denote the $Lr$-norm for any real $r≥1$⁠. The following lemma is a direct consequence of the covariance inequality of Ibragimov [12] and the inequality $2α≤β$⁠.

The weak consistency of the classifier (1.2) has been established by [21]. In this section we study the strong consistency of (1.2). The following theorem states the strong consistency under mild conditions.

Our aim in this section is to look at how the classifier (1.2) behaves on simulated samples by comparing it with the classical kernel rule. We use the R statistical programming environment to run a simulation study for $N=2$⁠. Let ${(X(i,j),Y(i,j))}$ be the field of interest and suppose that the simulated data are observed on the area $I(n,n)={(i,j)∈ℤ2:1≤i,j≤n}$⁠. Let

where $M={(2k,2l),1≤k,l≤10}$ is the set of non-observed sites which need to be classified. In this particular case, the vicinity of any missing site $(i,j)$ may be taken as in Figure 1.

It is important to note that the vicinity $ν(i,j)$ may be designed depending on the location of the missing site (see some typical examples in Figure 2) and that samples with larger size give more freedom to design vicinities.

Figure 2 shows some examples of vicinities that can be used when the missing sites are not completely surrounded by already labeled sites (located at the edges of $Sn$ for example).

We suppose that the simulated fields have the covariance function

We use the classifier (1.2) with $K(x)=∏i=18Ki(xi)$ for $x=(x1,…,x8)∈ℝ8$ where $Ki(xi)$ is the standard Gaussian density (Gaussian kernel). We suppose that ${X(i,j),1≤i,j≤n}$ are observations of a Gaussian mixture model:

with $μ0<μ1<μ2$ and $π1+π2+π3=1$⁠. In order to illustrate the fact that our method works for multi-class, the data set ${X(i,j),1≤i,j≤n}$ is partitioned in three clusters as follows:

For each $n=50,75,100$⁠, we generate $100$ samples on the region $I(n,n)$ with $μ0=5$⁠, $μ1=15$⁠, $μ2=25$⁠, $π0=π1=π2=1/3$ and $σ02=σ12=σ22=4$⁠. In each replication, we use the classifier (1.2), constructed on the basis of the training data observed on $J(n,n)$⁠, to re-predict the labels of sites in the test set $M$⁠. Figure 3 displays one replication for $n=50$⁠.

The optimal bandwidth $b^opt$ is obtained by minimizing the cross-validation criterion on a training sample and the misclassification error rate (⁠$ER$⁠) is evaluated based on the associated test sample. The average error rate (⁠$AER$⁠) is obtained by averaging the error rates associated with the corresponding $100$ test samples.

Table 1 shows that the estimated optimal bandwidth and the average error rate decrease when the training sample size increases. This means that the practical results in the simulation study are in line with the theoretical results. Now, let us compare the average error rate (⁠$AER$⁠) resulting from application of the proposed classifier with that resulting from application of the classical kernel rule.

The classical kernel rule is given, for any unlabeled site $j$ with $Xj=x$⁠, by

where $K˜:ℝd→ℝ+$ is a kernel on $ℝd$ (the Gaussian kernel is considered here), and $hn$ is a sequence of bandwidths. In order for the classical kernel classifier to be usable in our case, we have to adjust it slightly by taking the sum over $In−M$ instead of $In$⁠, $i.e.$⁠, for each $j∈M$ with $Xj=x$⁠,

From the theoretical point of view, this is justified by the fact that $g˜n$ has the same asymptotic behavior on $In$ as on $In−M$ since $M$ is bounded. In this classical kernel method, we consider knowing the features vector $Xj$ of each element $j$ of $M$ and we use $x$⁠, the value of $Xj$⁠, to predict its class while we needed only observations in nearby sites to predict the label of $j$ by the classifier (1.2). We apply the classical kernel classifier to re-classify the elements of $M$ using the same training samples generated above and taking into account all the replications for each size $n=50,75,100$⁠. Similar to what we have done in application of (1.2), the optimal bandwidth $h^opt$ is chosen by minimizing the cross-validation criterion on a training sample and the misclassification error rate (⁠$ER$⁠) is evaluated based on the associated test sample. Table 2 reports the average error rate (AER), obtained by averaging the error rates associated with the corresponding $100$ test samples.

By comparing Tables 1 and 2, we observe that the corresponding error values in the two tables begin to be close as $n$ increases. This supports the possibility of using the classifier (1.2) as an alternative to the classical kernel classifier when we have to classify sites with missing features.

A digital image is nothing than data numbers indicating variation of $red$⁠, $green$ and $blue$ (⁠$RGB$⁠) at a particular location on a grid of pixels. An $RGB$ color value is specified with: $rgb(red,green,blue)$⁠. Each parameter $(red,green,blue)$ defines the intensity of the color as an integer between $0$ and $255$⁠. For example, $rgb(0,0,255)$ is rendered as blue, because the blue parameter is set to its highest value $255$ and the others are set to $0$⁠. One can divide $RGB$ color values by $255$ in order to provide values in the interval $[0,1]$⁠. Let us have an image of Eiffel tower with $100$ missing pixels as in Figure 3.

We use the R package $jpeg$ to convert a $jpg$ image into $3$-d array of numbers. The package $jpeg$ offers the $readJPEG()$ function which can read raster graphics (consisting of “pixel matrices”) in $jpg$ format into $R$⁠. It returns either a single matrix with gray values in $[0,1]$ or $3$-d array with the $RGB$ values in $[0,1]$⁠, say $E$⁠. In our example of Figure 3, the dimensions of $E$ are $306×165×3$⁠. Thus, the elements of $E[,j]$ represent the intensities of the color $j$⁠, for $j=“red”,“green”or“blue”$⁠, at all pixels of the grid $I(306,165)$⁠. For example, the matrix $E[55:60,1:6,1]$ displays the intensities of $red$ in each pixel of the region:

Le $X(i,j)=(X(i,j)(1),X(i,j)(2),X(i,j)(3))$ where $X(i,j)(k)$ is the intensity of the color $k$ at the pixel $(i,j)$⁠. Since our purpose is to classify new sites with completely missing features, we set an arbitrary threshold of 0.4 and we define labels as follow:

The set of $100$ missing pixels is taken as a test set, say $M$⁠. We use the classifier (1.2) (see (1.7) for the binary version) to classify each element of $M$ based on its eight-neighbors. The optimal bandwidth is evaluated by minimizing the cross-validation criterion on the known sites where we get $b^opt≈0.72$⁠. The misclassification error rate (⁠$ER$⁠) is evaluated on $M$ where we obtain $ER=0.04$ which indicates that there are only four misclassified cases out of $100$ classified cases (see Figure 4).

Now let us use the support vector machine (S V M) classifier to re-classify the elements of $M$⁠. In this case we should suppose that the RGB value is known for each element of $M$⁠. For implementing support vector machine in R programming language, we use the package $e1071$⁠. According to this classifier, we get a misclassification error of $ER=0.11$ and this permits to conclude that our kernel classifier in this example proceeds well compared to the (SVM) procedure.

Without loss of generality, we prove the theorem in the binary case where $Yj$ takes values in ${0,1}$ since no additional argument is required to prove it in the multi-class case. However, the Bayes classifier (1.1) in the binary case is given by

and the classifier (1.2) is given by

Define

Consequently, the classifier (7.1) can be written as

By Theorem 2.3 in [7], the consistency will be proved if we show that

But

Hence, in order to prove (7.1), it suffices to show that

and

The proof of (7.3) is the same as in the i.i.d. case (see [7], pp. 156–157 ). So, it suffices to prove (7.4). To do that, we will employ the blocking technique used in [4]. Let $p=pn=[n^γ]$ for some $1/θ<γ<1/(2N)$ (where $[.]$ stands for the integer part). Without loss of generality, we suppose that there exists a positive integer $qk$ such that $nk=2pqk$ for each $k=1,…,N$⁠. Let

We define blocks as follow, for each $j∈Jq$⁠,

As a consequence, we have $In=⋃k=12N⋃j∈JqSj(k)$⁠, and for each $k=1,…,2N$⁠, $card(Sj(k))=pN$ and $dist(Sj(k),Sj′(k))≥p$ for any $j≠j′$⁠. Let $Γj(k)={i∈Sj(k):νi⊂Sn}$⁠, for each $k=1,…,2N$ and $j∈Jq$⁠. Hence, for a fixed $k$⁠, we have $dist(Γj(k),Γj′k)≥p$ for any $j≠j$⁠, $card(Γj(k))≤card(Sj(k))=pN$ and

Let ${(X(i)∗,Yi∗)}i∈In−Jn$ be a set of independent and identically distributed random vectors such that they are independent of ${(X(i),Yi)}i∈Jn$ and $(X(i)∗,Yi∗)$ is identically distributed with $(X(1),Y1)$⁠. In order to make sense to the blocking technique, we define random vectors as follow: for each $i∈In$⁠,

It is clear that ${(X(i),Yi),i∈Jn}={(X(i),Yi),i∈Jn}$ and ${(X(i),Yi),i∈Γj(k)}={(X(i),Yi),i∈Γj(k)}$⁠. Now, for a fixed $k$ and each $j∈Jq$⁠, let $Wj(k)={(X(i),Yi),i∈Sj(k)}$ be a vector whose components are ordered according to a given order on indices. Applying Lemma 3.2 together with the blocks decomposition introduced by [10] (see also [20]) on the family of vectors ${Wj(k),j∈Jq}$⁠, we can generate independent copies ${W˜j(k),j∈Jq}$ such that: they are mutually independent, and for each $j∈Jq$⁠, $W˜j(k)={(X˜(i),Y˜i),i∈Sj(k)}$ has the same distribution as $Wj(k)={(X˜(i),Y˜i),i∈Sj(k)}$⁠. Furthermore, by Lemma 3.5, we have $P(Wj(k)≠W˜j(k))≤ϕ(p−r˜)$ since $p≥r˜$ for $n^$ large enough. Thus, the two vectors $(X˜(i),Y˜(i))$ and $(X˜(i′),Y˜i′)$ are independent for each $i∈Sj(k)$ and $i′∈Sj′(k)$ with $j≠j′$⁠. Now, for each $i∈Jn$⁠, there exists $j∈Jq$ such that ${(X(i),Yi)≠(X˜(i),Y˜i)}⊆(Wj(k)≠W˜j(k))$⁠. Since $(X˜(i),Y˜i)=(X˜(i),Y˜i)$ for each $i∈Jn$⁠, denote $(X˜(i),Y˜i)=(X˜(i),Y˜i)$⁠, for each $i∈Jn$ (or $i∈Γj(k)$⁠). As a consequence

By (7.5), we can write