[R-sig-Geo] cokriging simulation

[Edzer Pebesma]

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