gstat now uses Lapack / failing cokriging
On 17/11/15 18:14, Bruin, Sytze de wrote:
Hi Edzer, thanks for your prompt reply. I tried reproducing your matrix r but got a different result, i.e. a valid covariance matrix. The cross-covariances are different. Using the same example as above, my code is as follows:
library(sp)
library(gstat)
# some data
x <- c(215, 330, 410, 470, 545)
y <- c(230, 310, 330, 340, 365)
fc <- c(0.211, 0.251, 0.281, 0.262, 0.242)
por <- c(0.438, 0.457, 0.419, 0.430, 0.468)
Allier <- data.frame(x, y, fc, por)
coordinates(Allier) = ~x+y
# gstat object for co-kriging using linear model of co-regionalization
g <- gstat(id=c("fc"), formula=fc~1, data=Allier,
model=vgm(0.00247, "Sph", 480, 0.00166))
g <- gstat(g, id="por", formula=por~1, data=Allier,
model=vgm(0.00239, "Sph", 480, 0.00118))
g <- gstat(g, id=c("fc", "por"),
model=vgm(0.00151, "Sph", 480, -0.00124))
dists <- spDists(Allier)
r11 <- variogramLine(g$model$fc, dist_vector = dists, covariance = T)
r12 <- variogramLine(g$model$fc.por, dist_vector = dists, covariance = T)
r22 <- variogramLine(g$model$por, dist_vector = dists, covariance = T)
r <- cbind(r11, r12)
r <- rbind(r, cbind(r12, r22))
r
[,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,] 0.0041300000 0.0014193874 0.0008959951 0.0005655821 0.0002240733 0.0002700000 0.0008677227
[2,] 0.0014193874 0.0041300000 0.0018397575 0.0013976206 0.0008790828 0.0008677227 0.0002700000
[3,] 0.0008959951 0.0018397575 0.0041300000 0.0020030001 0.0014238096 0.0005477541 0.0011247100
[4,] 0.0005655821 0.0013976206 0.0020030001 0.0041300000 0.0018652970 0.0003457607 0.0008544158
[5,] 0.0002240733 0.0008790828 0.0014238096 0.0018652970 0.0041300000 0.0001369841 0.0005374150
[6,] 0.0002700000 0.0008677227 0.0005477541 0.0003457607 0.0001369841 0.0035700000 0.0013734154
[7,] 0.0008677227 0.0002700000 0.0011247100 0.0008544158 0.0005374150 0.0013734154 0.0035700000
[8,] 0.0005477541 0.0011247100 0.0002700000 0.0012245061 0.0008704261 0.0008669750 0.0017801702
[9,] 0.0003457607 0.0008544158 0.0012245061 0.0002700000 0.0011403233 0.0005472637 0.0013523535
[10,] 0.0001369841 0.0005374150 0.0008704261 0.0011403233 0.0002700000 0.0002168159 0.0008506105
[,8] [,9] [,10]
[1,] 0.0005477541 0.0003457607 0.0001369841
[2,] 0.0011247100 0.0008544158 0.0005374150
[3,] 0.0002700000 0.0012245061 0.0008704261
[4,] 0.0012245061 0.0002700000 0.0011403233
[5,] 0.0008704261 0.0011403233 0.0002700000
[6,] 0.0008669750 0.0005472637 0.0002168159
[7,] 0.0017801702 0.0013523535 0.0008506105
[8,] 0.0035700000 0.0019381256 0.0013776943
[9,] 0.0019381256 0.0035700000 0.0018048825
[10,] 0.0013776943 0.0018048825 0.0035700000
eigen(r)$values
[1] 0.0124609506 0.0055040132 0.0046425761 0.0035980843 0.0031064016 0.0028989486 0.0028042439 [8] 0.0018107335 0.0010254131 0.0006486352
_______________________________________________ R-sig-Geo mailing list R-sig-Geo at r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Building upon your earlier example, predict(g, SpatialPoints(data.frame(x=450, y=350)), debug = 32) gives you the generalized covariance matrix that is used for the cokriging, which I looked at. gstat computes generalized covariances as C(0)-gamma(h) with C(0) = max(0, sum of the positive sill values), instead of the sill of all sill values. If one of the sill components is negative, this matters. I looked in the Ver Hoef & Cressie 1993 paper, but couldn't find out which one is right. Maybe Gerard can also take a look at it; the fix would be trivial.
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/20151117/d7d023d9/attachment.bin>