Thanks Gerard, that has now been fixed in
https://github.com/r-spatial/gstat/issues/106
The problem was probably introduced in 2016 when I swapped out the
shipped matrix library with R's native one. The two differed in how a
Choleski decomposition was returned .
On 27/04/2022 21:29, Heuvelink, Gerard wrote:
Dear List,
We wish to create simulations from a cokriging model. I used to be able to do
that in the past, but now I am struggling to get the correlations preserved. I
dimplified the problem to a simple case with cross-correlation 0.8 (see script
below), but simulations show no correlation at all. What do I do wrong? Thank
you for your help, Gerard
library(gstat)
set.seed(12345)
d <- data.frame(0,0,0,0)
names(d) <- c("x","y","lmCres","lmNres")
coordinates(d) = ~x+y
nugC <- 1; nugN <- 1; nugCN <- 0.8
psillC <- 5; psillN <- 5; psillCN <- 4
range <- 10
g_l <- gstat(id=c("lmCres"),formula=lmCres~1,data=d,beta=0,nmax=100,
model=vgm(psillC,"Sph",range=10,nugC))
g_l <- gstat(g_l,id="lmNres",formula=lmNres~1,data=d,beta=0,nmax=100,
model=vgm(psillN,"Sph",range=10,nugN))
g_l <-
gstat(g_l,id=c("lmCres","lmNres"),model=vgm(psillCN,"Sph",range=10,nugCN))
df <- data.frame(x=1,y=1)
coordinates(df) <- ~x+y
nsim <- 250
sim <- predict(g_l,df,nsim=nsim)
simdf <- as.data.frame(sim)
plot(as.numeric(simdf[1,3:(nsim+2)]),as.numeric(simdf[1,(nsim+3):((2*nsim)+2)]))
cor(as.numeric(simdf[1,3:(nsim+2)]),as.numeric(simdf[1,(nsim+3):((2*nsim)+2)]))
[[alternative HTML version deleted]]
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
--
Edzer Pebesma
Institute for Geoinformatics
Heisenbergstrasse 2, 48151 Muenster, Germany
Phone: +49 251 8333081
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
