require(geoR) set.seed(123) R <- grf(50, cov.pars=c(1, 0.25)) class(R) <- "geodata" plot(R) ## Y1 : mu1 = 50, sig1 = 5, phi=0.25 Y1 <- R Y1$data <- 50 + 5 * Y1$data plot(Y1) ## Y2 : mu1 = 2.3; sig1 = 0.3, phi=0.25 Y2 <- R Y2$data <- 2.3 + 0.3 * Y2$data plot(Y2) ## summary(Y1$data) summary(Y2$data) cor(Y1$data, Y2$data) plot(Y1$data, Y2$data) ## agora simando o termo independente para cada variavel ## tau1 = 2 ; tau2 = 0.08 Y1$data <- Y1$data + 2*rnorm(50) Y2$data <- Y2$data + 0.08*rnorm(50) summary(Y1$data) summary(Y2$data) cor(Y1$data, Y2$data) plot(Y1$data, Y2$data) ## ## Modelo com componentes espaciais individuais ## ## Y1 : mu1 = 50, sig01 = 5, phi=0.25, sig1 = 4, phi1 = 0.4 Y1 <- R Y1$data <- 50 + 5 * Y1$data + grf(grid=R$coords, cov.pars=c(16, 0.4))$data plot(Y1) ## Y2 : mu1 = 2.3; sig1 = 0.3, phi=0.25, sig2 = 0.2, phi2 = 0.15 Y2 <- R Y2$data <- 2.3 + 0.3 * Y2$data + grf(grid=R$coords, cov.pars=c(0.04, 0.15))$data plot(Y2) summary(Y1$data) summary(Y2$data) cor(Y1$data, Y2$data) plot(Y1$data, Y2$data) ## variogramas empiricos e teoricos plot(variog(Y1)) lines.variomodel(seq(0, 1.5, len=100), cov.model=c("exp","exp"), cov.pars=rbind(c(25, 0.25), c(16, 0.4)), nug=0) plot(variog(Y2)) lines.variomodel(seq(0, 1.5, len=100), cov.model=c("exp","exp"), cov.pars=rbind(c(0.09, 0.25), c(0.04, 0.15)), nug=0) ### ## simulando, estimando e predizendo ### 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] par(mfrow=c(1,2)) locs <- expand.grid((0:50)/30, (0:50)/30) pred1 <- predict(fit12, loc=locs) image(pred1) contour(pred1, add=T) pred2 <- predict(fit12, loc=locs, var=2) image(pred2) contour(pred2, add=T)