[R-sig-Geo] covariance matrix for kriging predictions

[Benedikt Gräler]

Dear Jui-Han,

I believe the following script does what you are looking after. The idea
is to start with the covariance matrix purely defined by the variogram
model (i.e. with the joint sill on the diagonal). Rescale this
covariance matrix first to correlations and then back to covariances
using the marginal variances obtained from the kriging predictor. I will
correct the krige0 function in gstat accordingly.

Let me know if this is what you have been looking for.

Best wishes,

 Ben

#--R-Script--#
# get some data
library(gstat)
library(sp)

data(meuse)
coordinates(meuse) = ~x+y
data(meuse.grid)
gridded(meuse.grid) = ~x+y

m <- vgm(0.59,"Sph", 897, 0.05)

# do the prediction to get the marginal kriging
# variance for each grid cell
pred <- krige(log(zinc)~1, meuse, meuse.grid, m)

# calculate the covariacne strucutre of the desired
# grid: 3103 x 3103 matrix
covMat <- variogramLine(m, dist_vector = spDists(meuse.grid,
meuse.grid), covariance = TRUE)

# rescale to correlation matrix
corMat <- cov2cor(covMat)

# scale back to the desired variances on the
# diagonal, just the same way backwards as cov2cor
# does it:
covMat <- corMat*matrix(sqrt(pred$var1.var) %x%
sqrt(pred$var1.var),3103,3103)
covMat[1:5,1:5]

# check for correlation
cov2cor(covMat[1:5,1:5])

# check prediction variance:
diag(covMat) - pred$var1.var

#--R-Script--#

On 12.12.2013 22:54, Edzer Pebesma wrote:
> 
> 
> On 12/12/2013 08:20 PM, Jui-Han Chang wrote:
>> Hi all,
>>
>> I am trying to get covariance matrix for kriging predictions. I used krige0
>> function in the gstat library with options fullCovariance=TRUE and
>> computeVar=TRUE. I got a covariance matrix but the values seem to be too
>> large. I converted the covariance matrix to correlation matrix and got
>> values larger than 1. Following are the codes to derive the matrices:
>>
>> library(gstat)
>> data(meuse)
>> coordinates(meuse) = ~x+y
>> data(meuse.grid)
>> gridded(meuse.grid) = ~x+y
>>
>> m <- vgm(.59, "Sph", 874, .04)
>> x <- krige0(log(zinc)~1, meuse, meuse.grid, model =
>> m,fullCovariance=TRUE,computeVar=TRUE)
>>
>> x$var[1:5,1:5]
>>
>>        [,1]      [,2]      [,3]      [,4]      [,5]
>> 1 0.3108421 0.2794621 0.2883862 0.3019975 0.2528104
>> 2 0.2794621 0.2414045 0.2546174 0.2727450 0.2074984
>> 3 0.2883862 0.2546174 0.2631935 0.2769452 0.2263174
>> 4 0.3019975 0.2727450 0.2769452 0.2863739 0.2499542
>> 5 0.2528104 0.2074984 0.2263174 0.2499542 0.1648699
>>
>> cov2cor(x$var[1:5,1:5])
>>
>>       [,1]     [,2]     [,3]     [,4]     [,5]
>> 1 1.000000 1.020188 1.008247 1.012201 1.116746
>> 2 1.020188 1.000000 1.010131 1.037331 1.040091
>> 3 1.008247 1.010131 1.000000 1.008764 1.086450
>> 4 1.012201 1.037331 1.008764 1.000000 1.150331
>> 5 1.116746 1.040091 1.086450 1.150331 1.000000
>>
>> I am not sure what I am getting here now. Is this covariance matrix for
>> kriging predictions? If not does anyone know how to get a covariance
>> matrix? Any help will be appreciated.
> 
> Thanks, this is clearly wrong. My feeling is that c0 in the computation
> of the kriging predictions covariance matrix needs be specified
> differently, not as a constant, however I haven't sorted out how, so far.
> 
> 
> 
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> 

-- 
Benedikt Gräler

ifgi - Institute for Geoinformatics
University of Muenster

http://ifgi.uni-muenster.de/graeler

Phone: +49 251 83-33082
Mail: ben.graeler at uni-muenster.de

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