[R-sig-Geo] Ordinary kriging variance of the prediction with error structure

[Edzer Pebesma]

Manuel,

gstat does kriging (exact interpolatoin) when specifying a nugget
variance, and smoothing when defining this nugget as an Err (measurement
error) term instead:

library(sp)
library(gstat)
demo(meuse, ask = FALSE, echo = FALSE)
krige(log(zinc)~1, meuse, meuse[1,],vgm(.6, "Sph", 900, .05))
krige(log(zinc)~1, meuse, meuse[1,],vgm(.6, "Sph", 900, Err = .05))

# From:

krige(log(zinc)~1, meuse[1:5,], meuse[1,],
   vgm(.6, "Sph", 900, .05), set = list(debug = 32))
krige(log(zinc)~1, meuse[1:5,], meuse[1,],
   vgm(.6, "Sph", 900, Err = .05), set = list(debug = 32))

You'll see that the difference is in the right-hand-side, the covariance
between observations and prediction location (C0).

Hth,


On 05/02/16 11:48, Antonio Manuel Moreno Ródenas wrote:
> Hello Ruben,
> 
> Thanks for your answer,
> I tried to check your reference, but still I'm having the same issue,
> 
> Cressie 1993, he says in 3.2.27:
> var(S|Y) = Sum(w_i)*V_io + m - c_me
> 
> Being:
> [w_i, vector of weights]
> [V_io, Semivariogram value from the sampled point to the predictied]
> [m, lagrange multiplier]
> [c_me, nugget (due to measurement error)]
> 
> This is equivalent to the notation I posted previously as C(h>0)= partial
> variance + nugget - V(h>0)
> [C, Covariance function]
> [V, Semivariogram function]
> 
> Substituting V_io and knowing that Sum(w_i)= 1, this gives:
> 
> Var(S| Y) = partial variance - Sum(w_i) * C_oi - m
> [C_oi, vector Covariance between predicted and sampled coordinates]
> 
> This is the expression that I'm working with, in which the nugget should
> not appear (unless x_o=x_i) isn't it?
> 
> I'm working with the covariance matrix as the nugget will only affect
> C(h=0) and not C(h>0) which is a property I need. And then the layout of
> the kriging system should be the same yet adding a nugget term in the
> diagonal of the covariance matrix.
> 
> Am I misunderstanding something?
> 
> Thanks and kind regards
> Antonio
> 
> 
> On 5 February 2016 at 10:56, rubenfcasal <rubenfcasal at gmail.com> wrote:
> 
>> Kriging the noiseless version of Y is not “the solution of the
>> (standard) kriging system with a nugget effect in the Covariance
>> structure” (the semivariances/variogram at lag 0 may not be zero).
>> See e.g. Cressie, 1993, p. 128 (for instance, eq. 3.2.27 shows the
>> correct expression of the kriging variance).
>>
>> Best regards, Ruben.
>>
>> El 04/02/2016 a las 15:56, Antonio Manuel Moreno Ródenas escribió:
>>> Dear r-sig-geo community,
>>>
>>> I would like to bring a conceptual question on the implementation of
>>> ordinary kriging in gstat.
>>>
>>> I'm trying to account for my measurement error in a OK scheme. I assume
>>> that my sampled vector Y(x_i) is a noisy realisation of S(x_i) (the real
>>> variable), thus:
>>> Y(x_i) = S(x_i) + e_i, where e_i is ~N(0,tau^2).
>>>
>>> If that error is assumed to follow a certain set of conditions (unbiased,
>>> uncorrelated between itself/the variable and tau=constant), this is
>>> analogous to the solution of the kriging system with a nugget effect in
>> the
>>> Covariance structure.
>>>
>>> I coded the kriging system and its solution. In order to assess if my
>>> implementation is correct I contrasted it with the krige function in
>> gstat.
>>> The predicted value at each point is the same, meaning that I got
>> correctly
>>> the weights of the system.
>>>
>>> However, I'm really confused when dealing with the variance in the
>>> prediction. Which should have this form:
>>>
>>> Var(S(x_o) | Y) = Var(S(x_o)) - w' * C_oi - mu
>>>
>>> w' weights vector
>>> C_oi vector Covariance between predicted and sampled coordinates
>>> mu lagrange multiplier
>>>
>>> If my objective would be to predict the signal of the variable S(x_o),
>> the
>>> term of Var(S(x_o)) will correspond to the partial variance, (*the sill
>>> without the nugget*). This is what I understood by following the notation
>>> of Model-based Geostatistics from Peter J.Diggle, where it is explicitly
>>> mentioned in (pag 137 (6.8)).
>>>
>>> However, I only get agreement in my comparison with the krige (gstat)
>>> variance results if I use the total sill as Var(S(x_o)) that is (partial
>>> variance + nugget). So my question is:
>>>
>>> I'am right by thinking that still Var(S(x_o)) should not include the
>> nugget?
>>> What is the outcome of krige in gstat when you consider a nugget? is it
>> the
>>> prediction of the signal? or is it the prediction of what you would
>> measure
>>> at that location?
>>>
>>> Kind regards,
>>>
>>> Antonio
>>>
>>>       [[alternative HTML version deleted]]
>>>
>>> _______________________________________________
>>> R-sig-Geo mailing list
>>> R-sig-Geo at r-project.org
>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>
>>
>>
>>         [[alternative HTML version deleted]]
>>
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> R-sig-Geo at r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
> 
> 	[[alternative HTML version deleted]]
> 
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> 

-- 
Edzer Pebesma
Institute for Geoinformatics  (ifgi),  University of Münster
Heisenbergstraße 2, 48149 Münster, Germany; +49 251 83 33081
Journal of Statistical Software:   http://www.jstatsoft.org/
Computers & Geosciences:   http://elsevier.com/locate/cageo/
Spatial Statistics Society http://www.spatialstatistics.info

-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 490 bytes
Desc: OpenPGP digital signature
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20160205/7aec5e65/attachment.bin>