# More Efficient Prediction for Ordinary Kriging to Solve a Problem in the Structure of Some Random Fields

Authors: Mohammad Mehdi Saber; Ramy Abdelhamid Aldallal
Journal: Complexity
 (2022)

---

## Abstract

Recently, some specific random fields have been defined based on multivariate distributions. This paper will show that almost all these random fields have a deficiency in spatial autocorrelation structure. The paper recommends a method for coping with this problem. Another application of these random fields is spatial data prediction, and the Kriging estimator is the most widely used method that does not require defining the mentioned random fields. Although it is an unbiased estimator with a minimum mean-squared error, it does not necessarily have a minimum mean-squared error in the class of all linear estimators. In this work, a biased estimator is introduced with less mean-squared error than the Kriging estimator under some conditions. Asymptotic behavior of its basic component will be investigated too.


---

## Body

## 1. Introduction

Random fields (RFs) as a statistical model can be applied in many real-life events, such as biological sequences, management of soil resources in agriculture and forestry, text and image processing, designing environmental monitoring networks, artificial intelligence, and road and tunnel planning.An RF is usually constructed based on a multivariate distribution’s finite-dimensional distributions. For example, the skew-normal (SN) RF based on multivariate SN distribution [1] has been defined by Kim and Mallick [2]. The closed skew-normal (CSN) RF based on multivariate CSN distribution [GF_2004] has been established by Allard and Naveau [3]. Similar work has been performed by Karimi et al [4], Hosseini et al. [5], and Karimi and Mohammadzadeh [6]. Zareifard and Jafari Khaledi [7] defined a second-order stationary RF named unified skew-normal (SUN) RF by using multivariate SUN distribution [8]. The generalized Skew-Normal (GSN) RF introduced by Mahmoudian [9] has also been constructed from multivariate GSN distribution as shown by Sahu et al. [10], while the generalized asymmetric Laplace (GAL) RF is defined by multivariate GAL distribution as shown by Kozubowski et al. [11]and was introduced by Saber et al. [12]. Some other RFs may be defined by multivariate skew-t distribution and its general form or multivariate extended skew-t (EST) distribution as shown by Arellano-Valle and Genton [13].For the spatial structure of the model, Wall [14] studied the spatial structure implied by the CAR and SAR models. According to Saber et al. [12] and Mahmoudian [15], almost all mentioned RFs in the previous paragraph are not well defined by the Kolmogorov existence theorem. Saber et al. [12] showed that these RFs do not consider proper spatial autocorrelation structure (PSAS). In the present article, we will show that almost all aforementioned skew RFs have this deficiency. In fact, they do not consider PSAS.A final application of any random field is spatial data prediction, which has been applied in some works of authors such as Basu and Reinsel [16], Saber and Nematollahi [17], and Saber [18]. Many researchers have studied the prediction of spatial data using mentioned unsuitable RFs. To solve the previously discussed problem, some authors such as Mojiri et al. [19] have considered a univariate distribution for detection of errors of the process, which results in a new problem in modeling spatial correlation. Some others such as Hosseini and Karimi [21] have used an approximate skew-normal RF. Although this work has some benefits for reaching a high percentage of PSAS, it has some fundamental problems, and it cannot solve the two mentioned problems.Kriging is a spatial prediction methodology and is the most widely used method of spatial data prediction that has increased during recent years. Because this method does not require defining the previously mentioned random fields that have problem in PSAS, Cressie [22] presented an unbiased linear estimator for Kriging and has a minimum mean-squared error (MSE), although it is not the best estimator in the larger class of all linear estimators. On the other hand, according to Moyeed and Papritz [23] linear Kriging is as good as any nonlinear method, especially for symmetric data. However, for skewed data, some of the nonlinear methods perform better in estimating prediction uncertainty. The focus of this article is on linear Kriging as we suggested a biased estimator for Kriging where the only interest is to minimize the MSE. Our recommended estimator has less MSE than the Cressie estimator under some conditions. Therefore, it is appropriate to use it for skewed data when we have a problem in defining a valid RF and for ordinary Kriging that does not have precision for nonlinear methods.The paper is organized as follows: In Section2, some requirements of the previously mentioned multivariate distributions are reviewed to show defined RFs that do not consider the PSAS property. After that, an estimator for the mean of the process and then a predictor with a minimum MSE in all the linear predictors for ordinary Kriging are presented in Section 3. Finally, Section 4 is devoted to investigating some asymptotic behavior for the recommended estimator.


## 2. A Problem in Some RFs

The first part of this section reviews the variance matrix of multivariate SN, CSN, EST, GSN, SUN, and GAL distributions.Consider ap-dimensional continuous random vector X=X1,…,Xp′that has an SN distribution, denoted by X∼SNpμ,Σ,λ. Then, we have(1)VarX=Σ−2π1+λTλΣ1/2λλTΣ1/2,and consider another p-dimensional continuous random vector X=X1,…,Xp′ but with a CSN distribution, denoted by X∼CSNp,qμ,Σ,Γ,ξ,Δ. Then, we have(2)VarX=Σ+ΣΓTΛΓΣ−ΣΓTΨΨTΓΣ,where Λ and Ψ are complicated matrices based on Φq0;ξ,Δ+ΓΣΓ′. The detail about these matrices is not essential for us, although they are present in the study conducted by González-Farías et al. [24].The variance matrix of a multivariate EST distribution is denoted byX∼ESTpμ,Σ,λ,ν,τ, which is given as follows:(3)VarX=α∗2ω−1Σω−1−δδT+σ∗2δδT.However, the variance matrix of a multivariate GSN distribution is denoted byX∼GSNpμ,Σ,λ, which is given as follows:(4)VarX=Σ+1−2πdiagλ2.Arellano-Valle and Azzalini [25] recently provided a very complicated and cumbersome computation of the variance matrix of multivariate SUN. Let X∼SUNd,mξ,Σ,Δ,τ,Γ. Then we have,(5)VarX=Σ+ωΔHΔTω.Finally, the variance matrix of the GAL distribution continuous p-dimensional random vectorX denoted by X∼GALpμ,Σ,q is given as follows:(6)VarX=qΣ+μμT.After that, we state the following remark.Remark 1.
LetZs1,…,Zsn be the observations from an RF Zs:s∈D⊆Rd at n locations s1,…,sn in the area D. It is assumed that CovZsi,Zsj=qsi−sj where h is the norm of the vector h and q is a real value function. In a matrix form, VarZ=C where Z=Zs1,…,Zsn and C=Ci,ji,j=1n with a notation Ci,j=qsi−sj. Some RFs have been defined so that they satisfy the condition VarZ=C, while in some other cases, the covariance matrix C has been used somehow in which VarZ≠C.
In continuation of this section, we investigate the condition whereVarZ=C is named as PSAS in some well-known RFs.
Gaussian RF:Z∼Nnμ,Σ and VarZ=Σ. Therefore, we choose Σ=C to have a spatial correlation between variables.
t RF: Z∼tnμ,Σ,dd≥3 and VarZ=d/d−2Σ . It seems to be an appropriate choice to reach the PSAS when Σ=d−2/dC.
SN RF:Z∼SNnμ,Σ,λ and arZ=Σ−2/π1+λTλΣ1/2λλTΣ1/2 where it is convenient to substitute Σ−2/π1+λTλΣ1/2λλTΣ1/2=C to get to the PSAS.
CSN RF:Z∼CSNn,qμ,Σ,Γ,ξ,Δ and arZ=Σ+ΣΓTΛΓΣ−ΣΓTΨΨTΓΣ . Choosing Σ+ΣΓTΛΓΣ−ΣΓTΨΨTΓΣ=C provides an adequate choice to reach the PSAS.
EST RF:Z∼ESTnμ,Σ,λ,ν,τ and arZ=α∗2ω−1Σω−1−δδT+σ∗2δδT. By substituting α∗2ω−1Σω−1−δδT+σ∗2δδT=C, we will obtain the PSAS.
SUN RF:Z∼SUNn,mξ,Σ,Δ,τ,Γ and arZ=Σ+ωΔHΔTω . Making Σ+ωΔHΔTω=C will fulfill the PSAS.
GSN RF:Z∼GSNnμ,Σ,λ1 and arZ=Σ+λ1−2/πI . It provides a suitable choice to choose Σ+λ1−2/πI=C and then consider λ as a nugget and Σ+λ1−2/πI=C to accomplish the PSAS.
GAL RF:Z∼GALnμ,Σ,q and arZ=qΣ+μμT . After finding the PSAS, we will choose qΣ+μ211T=C, which results as μ=Xβ and Σ=1/qC−XββTXT where X is the matrix of covariates and β is the regression coefficient.
As previously shown, some of the RFs satisfy the PSAS, while others do not. In the following remark, we denote this matter concisely.Remark 2.
The Gaussian RF,t RF, GSN RF, and GAL RF satisfy the PSAS. Other RFs containing the SN RF, CSN RF, EST RF, and SUN RF do not satisfy the PSAS. In some of the previous works that applied a second group of RFs, it was supposed that Σ=C, although it does not lead to VarZ=C for achieving the PSAS.
Imposing the equalityVarZ=C, without any other constraints, may result in a nonpositive definite for the matrix Σ. There is a method from which we can add a small value to the diagonal of Σ, but it does not work in every case. So, there is still a problem in getting the PSAS and positive definiteness for the matrix. Σ


## 3. Linear Kriging as an Alternative for the Nonlinear Predictors Conducted by RFs

Let RFZs:s∈D⊆Rd satisfies the following model:(7)Zs=μ+Ys+εs,where μ is the mean value, Y (·) is the zero-mean intermediate-scale variation, and a white noise error denoted by ε (·) with a variance σε2.Kriging belongs to the family of linear least square estimation algorithms. Given the observationsZ = (Z (s1),…, Z(sn)), a common task is to predict Zs0 at an unobserved location s0 and to calculate the prediction error variance at each such location. The Kriging estimator is a linear predictor for Z^ (s0), which is as follows:(8)Z^so=∑i=1nliZsi,where lis (the Kriging weights) are some scalars that may be found under two conditions:(a)
EZ^s0−Zs0=0 (unbiased condition)(b)
ArgminMSEZ^s0=EZ^s0−Zs02Simple Kriging (known mean) and ordinary Kriging (unknown but constant mean) are two well-known divisions of Kriging based on the mean value specification. In fact, in simple Krigingμ is assumed to be known, but it is unknown in the case of ordinary Kriging, which occurs more in the application. As in the previous section, let C=covZsi,Zsji,j denotes the covariance matrix of the observations and let γ=covZs0,Zs1,…,covZs0,Zsnˊ denotes the vector containing the covariances between observed and unobserved location s0.According to Cressie [22], the simple Kriging estimator and its error variance are given, respectively, by(9)Z^simples0=γˊC−1Z,and(10)VarZ^simples0−Zs0=σz2s0−γˊC−1γ.Also, the ordinary Kriging and its error variance are given, respectively, by(11)Z^s0=γ+11−1ˊC−1γ1ˊC−11ˊC−1Z,and(12)σr2s0=VZ^s0−Zs0=σz2s0−γˊC−1γ+1−1ˊC−1γ21ˊC−11,where 1ˊ=1,…,1 is the row vector of size n in which its components are 1s.


### 3.1. New Estimator

The Kriging estimator for ordinary Kriging is the best estimator regarding unbiasedness and the minimum MSE. However, a biased estimator with a less MSE is preferable in some circumstances. An estimator with this condition exists, and we denote our motivation for recommending this estimator in the following theorem.Theorem 1.
An unbiased linear estimator for the mean of the process given in (7) exists and has a minimum MSE among all unbiased linear estimators for the mean. This estimator is in the form of(13)μ^=1ˊC−1Z1ˊC−11,and its variance is given by(14)Varμ^=11ˊC−11.Proof.
A linear estimator forμ is in the form of ∑i=1naiZsi=aˊZ. Unbiasedness implies that aˊ1=∑i=1nai=1. Hence, under this condition, VaraˊZ=aˊCa must be minimized. By using the Lagrange equation Γ=aˊCa+λaˊ1−1 and by vector differentiating, one can solve this problem to reach a^=C−11/1ˊC−11, and its results are shown in (13). Now, we haveVar1ˊC−1/1ˊC−11Z=1ˊC−1CC−11/1ˊC−112=1ˊC−11/1ˊC−112=1/1ˊC−11, which completes the proof.
Now, a simple Kriging method can be applied instead of an ordinary Kriging method. By substituting (13) in (7), we can reach Zs−μ^=Ys+εs and Zs−1ˊC−1Z/1ˊC−11=Ys+εs. From (9), we can conclude that(15)Z^news0=γˊC−1Z−1ˊC−1Z1ˊC−11=γ−11ˊC−11ˊC−1Z,which is the new estimator for interpolation. Z^news0 is a special case of Z^s0 when γˊC−11=2.Corollary 1.
The biasedness and prediction error variance of the estimatorZ^new is given by(16)EZ^news0−Zs0=μγˊC−11−2,VarZ^news0−Zs0=σz2s0−γˊC−1γ+11ˊC−11.Proof.
The proof comes from Theorem 1, and a direct computation is presented in the following equations:(17)EZ^news0−Zs0=γ−11ˊC−11ˊC−1μ1−μ=μγˊC−11−2,VarZ^news0−Zs0=VarZ^news0+VarZs0−2covZ^news0,Zs0=γ−11ˊC−11ˊC−1γ−11ˊC−11+σz2s0−2γ−11ˊC−11ˊC−1γ=σz2s0−γˊC−1γ+11ˊC−11.
Because of this fact, we deduce thatZ^news0 is a biased estimator. Now, to compare with the ordinary Kriging estimator given by (11), variance is not a useful criterion. The MSE is better when solving this problem and not only in this case but also while comparing between estimators even though at least one of them is biased. Finally, Theorem 2 shows that Z^news0 is better than the ordinary Kriging estimator (11) regarding the MSE under some conditions.Theorem 2.
Under one of the following two conditions (i) and (ii),MSE^Z^news0 will be less than MSE^Z^s0.(i)
1ˊC−1Z2>1ˊC−11 and γˊC−11 are between 2 and 21ˊC−1Z2/1ˊC−1Z2−1ˊC−11(ii)
1ˊC−1Z2<1ˊC−11 and γˊC−11>2 or γˊC−11<21ˊC−1Z2/1ˊC−1Z2−1ˊC−11)Proof.
From Corollary1, we have MSEZ^news0=μγˊC−11−22+σz2s0−γˊC−1γ+1/1ˊC−11 . On the other hand, Equation (12) results in MSEZ^s0=σz2s0−γˊC−1γ+1−γˊC−112/1ˊC−11. Then, by defining c=γˊC−11 and a=1ˊC−11, one can show that d=MSEZ^news0−MSEZ^s0=2c−c2/a+μ22−c2
From the positive definiteness ofC, it is found that a>0. By multiplying the above equality by a, we get ad=2c−c2+aμ22−c2. Replacing μ by μ^ leads to(18)ad^=c2e−a+c2a−4e+4e,where e=1ˊC−1Z2. The quadratic form (18) has two roots c1=2 and c2=2e/e−a with respect to c. Under condition (i) or (ii), Equation (18) will be negative; hence, d^≤0 under condition (i) or (ii) as well. This completes the proof.


## 3.1. New Estimator

The Kriging estimator for ordinary Kriging is the best estimator regarding unbiasedness and the minimum MSE. However, a biased estimator with a less MSE is preferable in some circumstances. An estimator with this condition exists, and we denote our motivation for recommending this estimator in the following theorem.Theorem 1.
An unbiased linear estimator for the mean of the process given in (7) exists and has a minimum MSE among all unbiased linear estimators for the mean. This estimator is in the form of(13)μ^=1ˊC−1Z1ˊC−11,and its variance is given by(14)Varμ^=11ˊC−11.Proof.
A linear estimator forμ is in the form of ∑i=1naiZsi=aˊZ. Unbiasedness implies that aˊ1=∑i=1nai=1. Hence, under this condition, VaraˊZ=aˊCa must be minimized. By using the Lagrange equation Γ=aˊCa+λaˊ1−1 and by vector differentiating, one can solve this problem to reach a^=C−11/1ˊC−11, and its results are shown in (13). Now, we haveVar1ˊC−1/1ˊC−11Z=1ˊC−1CC−11/1ˊC−112=1ˊC−11/1ˊC−112=1/1ˊC−11, which completes the proof.
Now, a simple Kriging method can be applied instead of an ordinary Kriging method. By substituting (13) in (7), we can reach Zs−μ^=Ys+εs and Zs−1ˊC−1Z/1ˊC−11=Ys+εs. From (9), we can conclude that(15)Z^news0=γˊC−1Z−1ˊC−1Z1ˊC−11=γ−11ˊC−11ˊC−1Z,which is the new estimator for interpolation. Z^news0 is a special case of Z^s0 when γˊC−11=2.Corollary 1.
The biasedness and prediction error variance of the estimatorZ^new is given by(16)EZ^news0−Zs0=μγˊC−11−2,VarZ^news0−Zs0=σz2s0−γˊC−1γ+11ˊC−11.Proof.
The proof comes from Theorem 1, and a direct computation is presented in the following equations:(17)EZ^news0−Zs0=γ−11ˊC−11ˊC−1μ1−μ=μγˊC−11−2,VarZ^news0−Zs0=VarZ^news0+VarZs0−2covZ^news0,Zs0=γ−11ˊC−11ˊC−1γ−11ˊC−11+σz2s0−2γ−11ˊC−11ˊC−1γ=σz2s0−γˊC−1γ+11ˊC−11.
Because of this fact, we deduce thatZ^news0 is a biased estimator. Now, to compare with the ordinary Kriging estimator given by (11), variance is not a useful criterion. The MSE is better when solving this problem and not only in this case but also while comparing between estimators even though at least one of them is biased. Finally, Theorem 2 shows that Z^news0 is better than the ordinary Kriging estimator (11) regarding the MSE under some conditions.Theorem 2.
Under one of the following two conditions (i) and (ii),MSE^Z^news0 will be less than MSE^Z^s0.(i)
1ˊC−1Z2>1ˊC−11 and γˊC−11 are between 2 and 21ˊC−1Z2/1ˊC−1Z2−1ˊC−11(ii)
1ˊC−1Z2<1ˊC−11 and γˊC−11>2 or γˊC−11<21ˊC−1Z2/1ˊC−1Z2−1ˊC−11)Proof.
From Corollary1, we have MSEZ^news0=μγˊC−11−22+σz2s0−γˊC−1γ+1/1ˊC−11 . On the other hand, Equation (12) results in MSEZ^s0=σz2s0−γˊC−1γ+1−γˊC−112/1ˊC−11. Then, by defining c=γˊC−11 and a=1ˊC−11, one can show that d=MSEZ^news0−MSEZ^s0=2c−c2/a+μ22−c2
From the positive definiteness ofC, it is found that a>0. By multiplying the above equality by a, we get ad=2c−c2+aμ22−c2. Replacing μ by μ^ leads to(18)ad^=c2e−a+c2a−4e+4e,where e=1ˊC−1Z2. The quadratic form (18) has two roots c1=2 and c2=2e/e−a with respect to c. Under condition (i) or (ii), Equation (18) will be negative; hence, d^≤0 under condition (i) or (ii) as well. This completes the proof.


## 4. Asymptotic Behavior of the Mean Estimator

Since the estimator ofμ (13) in Theorem 1 has an essential role, its asymptotic behavior should be interesting for this study. First, we denote the following lemma, which is required for the continuation of this section.Lemma 1.
For two sequences of unbiased estimatorsθ^n1 and θ^n2 we can say that θ^n2⟶L2θ, θ^n1⟶L2θ , if Varθ^n2≥Varθ^n1. Here, ⟶L2 states the mean square convergence.Theorem 3.
Let a random processZ (s) exists on a lattice Λ in Rk with coordinates, s1,…,sn, VarZsi≤c<∞ for all i=1,…,n and covZsi,Zsj⟶0 as si−sj tends to infinity, then Z¯n⟶L2μ and μ^⟶L2μ.Proof.
Without losing generality, assumeZi denotes Zsi. First, we know that EZ¯n=1/n∑i=1nEZi=μ. Now, we can compute the variance of Z¯n(19)VarZ¯n=1n2∑i=1n∑j=1nVarZiVarZjρZi,Zj≤cn2∑i=1n∑j=1nρZi,Zj=cn2∑i=1n∑j=1nρZi,ZjIsi−sj≥Nε+cn2∑i=1n∑j=1nρZi,ZjIsi−sj≥Nε.
The last term hasn2 components. The assumption which ρZsi,Zsj⟶0 , as si−sj tends to infinity gives that for all ε>0, there exists an integer Nε enough large such that ρZsi,Zsj≤ε for all Zsi and Zsj and that si−sj≥Nε. For a fixed si, the number of sj so that si−sj≤Nε is less than 2Nεk. Therefore, the total number of points with condition si−sj≤Nε is at last n2Nεk which leads to VarZ¯n≤c/n2n2Nεk+n2−n2Nεkε=c1/n2Nεk+1−1/n2Nεkε
which clearly provesZ¯n⟶L2μ. By Theorem 1, μ^ is unbiased for μ and Varμ^<VarZ¯n. Finally, by Lemma 1 we get that μ^⟶L2μ.Corollary 2.
It is an immediate consequence of Theorem4 that under its conditions 1ˊC−11 tends to infinity as n tends to infinity.Theorem 4.
Under conditions of Theorem4 and at least one of the followings (i) or (ii) and Z^ns0⟶L2μ.(i)
γˊC−1γ⟶0 as n tends to infinity.(ii)
1ˊC−11γˊC−1γ−γˊC−112 is bounded.Proof.
Regarding Lemma1 and the unbiasedness of Z^ns0, we see that varZ^ns0⟶0 if it satisfies condition (i) or (ii). Since μn^ is the linear unbiased estimator for μ with the minimum variance and Z^ns0 is the linear unbiased estimator for μ, we obtain that(20)0≤varμn^−varZ^ns0=kn.
Corollary1 and equation (20) result in varZ^ns0⟶0 if and only if kn⟶0. On the other hand, a direct computation shows that kn=1ˊC−11γˊC−1γ−γˊC−112/1ˊC−11 and under condition (ii), it obviously tends to 0. Another form of kn is kn=γˊC−1γ−γˊC−112/1ˊC−11, with the fact that γˊC−112/1ˊC−11>0 and if (i) satisfies, it will tend to 0.


## 5. Conclusion

In this paper, we showed some previously defined RFs that have been extensively applied to modeling time series, spatial, and spatiotemporal data. However, they have a few limitations that they cannot consider the PSAS, which has a crucial and basic role in modeling the previously mentioned data. For that, we presented some alternative RFs without the mentioned problem, and having a more general solution is still needed. The paper recommended a method for coping with this problem by finding a new predictor in ordinary Kriging. The customary predictor for ordinary Kriging has a minimum MSE in the class of all linear unbiased predictors. However, it does not necessarily have a minimum MSE in the class of all linear predictors. Therefore, we obtained a biased predictor with a less MSE than the Kriging predictor. Asymptotic behavior of this predictor alongside the past Kriging predictor was provided.We studied six well-known multivariate distributions in this article. Some other multivariate distributions studied by Kotz et al. [26] such as the truncated multivariate normal distribution, truncated multivariate t distribution, Linnik’s distribution, gamma distribution, and logistic distribution can be surveyed under this paper’s point of view.Despite providing some solutions for the problem of lacking PSAS, it seems that the best method for modeling skew spatial data is defining a consistent skew RF. This matter may be performed in the future.



---

*Source: 9712576-2022-07-06.xml*