require(geoR) ## ## Testing code: Inference in the ## Bivariate Gaussian Common Component Model ## ## 1. Simulating Y1 and Y2 with some (but not all) locations in common ## ## 1.1. Coordinates set.seed(30) all.coords <- round(cbind(runif(200), runif(200)), dig=3) set.seed(111); ind1 <- sample(1:200, 120) coords1 <- all.coords[ind1,] set.seed(222); ind2 <- sample(1:200, 80) coords2 <- all.coords[ind2,] ## ## Model parameters: ## m1 = 10, mu2 = 50 ## sigma01 = 2, sigma1 = 1.5, sigma02 = 8, sigma2 = 6 ## phi0 = 0.2, phi1=0.15, phi2=0.25 mu1 <- 10; mu2 <- 50 sigma01 <- 2; sigma1 <- 1.5; sigma02 <- 8; sigma2 <- 6 phi0 <- 0.2; phi1 <- 0.15; phi2 <- 0.25 ## ## Model parameters (reparametrised) ## sigma = 2, nu1 = 0.75, eta = 4, sigma2 = 3 sigma <- sigma01 nu1 <- sigma1/sigma eta <- sigma02/sigma; nu2 <- sigma2/sigma ## ## Simulating model components ## S0 <- grf(grid=all.coords, cov.pars=c(sigma, phi0))$data S1 <- grf(grid=coords1, cov.pars=c(sigma1, phi1))$data S2 <- grf(grid=coords2, cov.pars=c(sigma2, phi2))$data ## ## Y1 and Y2 data ## y1 <- as.geodata(cbind(coords1,10+S0[ind1]+S1)) y2 <- as.geodata(cbind(coords2,50+S0[ind2]+S2)) ## ## maximum likelihood estimation ## fit12 <- likfitBGCCM(y1, y2, ini.s=c(2,1.5,8,6), ini.phi=c(.2,.15,.25), control=list(trace=T)) fit12 fit12a <- likfitBGCCM(y1, y2, ini.s=c(2,1.5,8,6), ini.phi=c(.2,.15,.25), method="BFGS") fit12a fit12b <- likfitBGCCM(y1, y2, ini.s=c(2,1.5,8,6), ini.phi=c(.2,.15,.25), method="L-BFGS-B", control=list(trace=T)) fit12b fit12c <- likfitBGCCM(y1, y2, ini.s=c(2,1.5,8,6), ini.phi=c(.2,.15,.25), fc.min="nlminb") fit12c ## ## Prediction ## locs <- cbind(c(0.2, 0.5, 0.7, 0.2), c(0.3, 0.5, 0.2, 0.8)) pred1 <- predict(fit12, loc=locs) pred1[1:2] pred2 <- predict(fit12, loc=locs, var=2) pred2[1:2]