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]]