#################################################### # An?lise para dados Parcialmente Co-locados # #################################################### rm(list=ls(all=TRUE)) require(geoR) LP <- list(nu1 = 0.9 , nu2 = 0.6, a = 1/4 , rho = 0.85 , sig1 = sqrt(0.7), sig2 = sqrt(0.5) , tau1 = sqrt(0.2), tau2 = sqrt(0.1) ) set.seed(29) gs <- cbind(runif(340, 0, 25), runif(340, 0, 25))#todas coordenadas gs dM <- dist(as.matrix(gs), upper=T, diag=T) dM lower2matrix <- function(x){ n <- (1+sqrt(1+8*length(x)))/2 M <- diag(n) M[lower.tri(M)] <- x M <- t(M) M[lower.tri(M)] <- x M } montaSigma <- function(parList, distPares){ evalq({ n <- (1+sqrt(1+8*length(distPares)))/2 ## condi??o/restri??o if(abs(rho) > sqrt(nu1 * nu2)/((nu1+nu2)/2)) stop("restri??es nos par?metros n?o satisfeitas") c11 <- cov.spatial(distPares,cov.model="matern",cov.pars=c(1, 1/a),kappa=nu1) c11 <- sig1^2 * lower2matrix(c11) diag(c11) <- diag(c11) + tau1^2 c22 <- cov.spatial(distPares,cov.model="matern",cov.pars=c(1, 1/a),kappa=nu2) c22 <- sig2^2 * lower2matrix(c22) diag(c22) <- diag(c22) + tau2^2 c12 <- cov.spatial(distPares,cov.model="matern",cov.pars=c(1, 1/a),kappa=(nu1+nu2)/2) c12 <- rho * sig1 * sig2 * lower2matrix(c12) S <- rbind(cbind(c11, c12), cbind(c12, c22)) }, env=parList) } dV <- as.vector(dist(gs)) system.time(SIGMA <- montaSigma(LP, distPares=dV)) cholSIGMA <- chol(SIGMA) var<- drop(crossprod(cholSIGMA, rnorm(dim(cholSIGMA)[1]))) var length(var) var <- var + c(rep(2, 340), rep(10, 340)) plot(variog(coords=gs, data=var[1:340]), ty="b") lines(variog(coords=gs, data=var[341:680]), ty="b", col=2) var1<-var[1:(length(var)/2)] var1 var2<-var[((length(var)/2)+1):length(var)] var2 ##################################################################### ####Para predi??o: 100pontos coords.pred<-cbind(gs[c(241:340)],gs[c(241:340),2]) var1.pred<-var1[c(241:340)] var2.pred<-var2[c(241:340)] ##################################################################### ####Utilizar na análise: var1<-var1[c(1:240)] var2<-var2[c(1:240)] coords<-cbind(gs[c(1:240)],gs[c(1:240),2]) DADOS<-cbind(coords,var1,var2) ###2)Dados parcialmente co-locados dados2<-cbind(DADOS[c(181:240),1],DADOS[c(181:240),2], DADOS[c(181:240),3],DADOS[c(181:240),4]) #60pontos co-locados x1<-dados2[,1] x2<-dados2[,2] var1<-dados2[,3] var2<-dados2[,4] #Var1: +60pontos b1<-cbind(DADOS[c(61:120),1],DADOS[c(61:120),2], DADOS[c(61:120),3],c(rep(NA,60))) #Var2: +60pontos b2<-cbind(DADOS[c(121:180),1],DADOS[c(121:180),2], c(rep(NA,60)),DADOS[c(121:180),4]) #CONJUNTO: x1<-c(x1,b1[,1],b2[,1]) x2<-c(x2,b1[,2],b2[,2]) var1.<-c(var1,b1[,3],b2[,3]) var2.<-c(var2,b1[,4],b2[,4]) DADOS2<-cbind(x1,x2,var1.,var2.) coords2<-cbind(DADOS2[,1],DADOS2[,2]) ################## BGCCM ####################### banco<-DADOS2 var1.BGCCM<-as.geodata(banco,coords.col=1:2,data.col=3) var2.BGCCM<-as.geodata(banco,coords.col=1:2,data.col=4) points(var1.BGCCM) plot(var1.BGCCM,low=T) points(var2.BGCCM) plot(var2.BGCCM,low=T) bgccm1<-likfitBGCCM(var1.BGCCM, var2.BGCCM, ini.sigmasq=c(8,4,5,2), ini.phi=c(6,3,4),cov0.model="matern", cov1.model="matern", cov2.model="matern", kappa0=0.5, kappa1=0.5, kappa2=0.5, control=list(trace=T)) bgccm1 dM. <- dist(as.matrix(coords2), upper=T, diag=T) dM. #da erro varcov<-varcovBGCCM(dM., cov0.pars=c(2,6), cov1.pars=c(3,7), cov2.pars=c(4,9), cov0.model = "matern", cov1.model = "matern", cov2.model = "matern", kappa0 = 0.4, kappa1 = 0.4, kappa2 = 0.6, scaled = TRUE, inv = FALSE, det = FALSE) varcov gr <- expand.grid(seq(0,1, by=0.05), seq(0,1, by=0.02)) #grid para predição predição1<-predict(object=bgccm1,locations=gr, variable.to.predict = 2) predição1 image(predição1, col=terrain.colors(15)) image(predição1, col=gray(seq(.8, 0, leng=15))) ####Predição dos valores retirados! pred1<-predict(object=bgccm1, locations=coords.pred,var=1) pred1 pred2<-predict(object=bgccm1, locations=coords.pred,var=2) pred2 image(pred1,col=terrain.colors(15)) #################### Gstat ##################### require(gstat) gs1<-DADOS2[,1] gs2<-DADOS2[,2] var1.gstat<-DADOS2[,3] var2.gstat<-DADOS2[,4] base<-as.data.frame(cbind(gs1,gs2,var1.gstat,var2.gstat)) base base.pred<-as.data.frame(cbind(coords.pred,var1.pred,var2.pred)) coordinates(base) <- ~gs1+gs2 #Para x1# variog_x1<-variogram(var1.gstat~1,data=base,width=0.9,cutoff=12) #variog_x1$gamma->quantifica??o da vari?ncia em uma dist?ncia h plot(variog_x1) mod1_x1<-vgm(psill=0.95,'Exp',range=2 ,nugget=0.1) fit.x1<-fit.variogram(variog_x1,mod1_x1,fit.sills = FALSE, fit.ranges = FALSE) plot(variog_x1,model=fit.x1) #Para x2# variog_x2<-variogram(var2.gstat~1,data=base,width=1,,cutoff=11) plot(variog_x2) mod1_x2<-vgm(psill=0.5,'Exp',range=3 ,nugget=0.2) fit.x2<-fit.variogram(variog_x2,mod1_x2,fit.sills = FALSE, fit.ranges = FALSE) plot(variog_x2,model=fit.x2) #Grid para Krigagem# x <- rep(seq(0,25,1),25) y <- rep(0:25,each=25,times=1) xy <- cbind(x,y) S <- SpatialPoints(xy) class(S) gridded(S) <- TRUE gridded(S) x1.krg<-krige(var1.gstat~1,base,S,fit.x1) image(x1.krg) summary(x1.krg) x2.krg<-krige(var2.gstat~1,base,S,fit.x2) image(x2.krg) summary(x2.krg) g1 <- gstat(NULL,id="x1",form=var1.gstat~1, data=base) g1 <- gstat(g1, id ="x2", form = var2.gstat~1, data=base) v.cross <- variogram(g1,width=0.7,cutoff=12) #variograma cruzado plot(v.cross) lmc<-fit.lmc(v=v.cross, g=g1, model=mod1_x1) lmc #Modelo Linear de Corregionalização plot(variogram(lmc),model=lmc$model) ####Predição dos valores retirados! pred.gstat1<-predict.gstat(lmc,base.pred) #não pred.gstat<-predict.gstat(lmc,base) pred.gstat ####################### sbBayes ################## require(spBayes) require(MASS) require(MBA) q<-2 #número de variáveis g<-180 #número de coordnadas m<-60 #número de nós nltr <- q*(q-1)/2+q var1.bayes<-var1 var2.bayes<-var2 n.samples<-3000 K.starting<-c(0.7,0.8,0.5) #matriz espacial de covariância cruzada Psi.starting<-c(0.2,0.15,0.1) #matriz não espacial de covariância cruzada m.1 <- spMvLM(list(var1.bayes~1,var2.bayes~1), coords=coords, starting=list("beta"=rep(1,q), "phi"=rep(3/0.5,q), "A"=K.starting, "L"=Psi.starting), sp.tuning=list("phi"=rep(0.1,q), "A"=rep(0.1,nltr), "L"=rep(0.1,nltr)), priors=list("phi.Unif"=rep(c(3/1,3/0.1),q), "K.IW"=list(q+1, diag(1,q)), "Psi.IW"=list(q+1, diag(1,q))), modified.pp=TRUE, cov.model="exponential", n.samples=n.samples, sub.samples=c(1500,n.samples,5), verbose=TRUE, n.report=100) m.1$p.samples[,paste("phi_",1:q,sep="")] <- 3/m.1$p.samples[,paste("phi_",1:q,sep="")] print(summary(mcmc(m.1$p.samples))) plot(m.1$p.samples) var.hat <- rowMeans(m.1$sp.effects) var.hat.1 <- var.hat[1:(length(var.hat)/2)] var.hat.2 <- var.hat[122:length(var.hat)] par(mfrow=c(2,2)) surf <- mba.surf(cbind(coords,var1.bayes), no.X=100, no.Y=100, extend=T)$xyz.est image(surf, main="Observed 1"); contour(surf, add=TRUE) surf <- mba.surf(cbind(coords,var.hat.1), no.X=100, no.Y=100, extend=T)$xyz.est image(surf, main="Fitted 1"); contour(surf, add=TRUE) points(m.1$knot.coords, pch=19, cex=1) surf <- mba.surf(cbind(coords,var2.bayes), no.X=100, no.Y=100, extend=T)$xyz.est image(surf, main="Observed 2"); contour(surf, add=TRUE) surf <- mba.surf(cbind(coords,var.hat.2), no.X=100, no.Y=100, extend=T)$xyz.est image(surf, main="Fitted 2"); contour(surf, add=TRUE) points(m.1$knot.coords, pch=19, cex=1) #################### ARTIGO ##################### geral<-function(pars1,pars2, printpars = FALSE, ...){ pars<-as.list(c(pars1,pars2)) rodaPARS <<- as.vector(pars) TODOS <- c("mu1","mu2", "nu1","nu2", "sig1","sig2","tau1","tau2","a","rho") if(any(is.na(match(TODOS, names(pars))))) stop("tá faltando parametro!!!") if(pars$nu1 <=0 || pars$nu2 <=0) return(-.Machine$double.xmax^0.5) # if(pars$sig1 <=0 || pars$sig2 <=0) return(-.Machine$double.xmax^0.5) ll<-function(pars, dados, logpars = FALSE, conc = FALSE) { # parsList recebe desvios (e n?o vari?ncias!) padr?o com op??o para log(desvio padrao) # dados : list(d, X, n1, n2, n)1e-10 if(!all(c("mu1","mu2","sig1") %in% names(pars1))) conc <- FALSE if(conc){ if(logpars) pars$sig2 <- exp(pars$sig2) pars$tau1 <- exp(pars$tau1) pars$tau2 <- exp(pars$tau2) pars$sig1 <- 1 } else{ if(logpars) pars$sig1 <- exp(pars$sig1) pars$sig2 <- exp(pars$sig2) pars$tau1 <- exp(pars$tau1) pars$tau2 <- exp(pars$tau2) } if(with(pars, abs(rho) > sqrt(nu1 * nu2)/((nu1+nu2)/2))) return(-.Machine$double.xmax^0.5) SIGMA <-montaSigma(pars, dados$d) tchSIGMA <- t(chol(SIGMA)) if(conc){ ## matrix (supondo n1 = n2) !!!!!! ONE <- kronecker(diag(1, 2), rep(1, dados$n1)) ## se n1 != n2: ## ONE <- cbind(c(rep(1, dados$n1), rep(0, dados$n2)), c(rep(0, dados$n1), rep(1, dados$n2))) ONEchSIGMA <- forwardsolve(tchSIGMA, ONE) mus <- solve(crossprod(ONEchSIGMA), crossprod(ONEchSIGMA, forwardsolve(tchSIGMA, dados$X))) } else mus <- c(pars$mu1,pars$mu2) #pars[1:2] Mu <- c(rep(mus[1],dados$n1), rep(mus[2],dados$n2)) res <- dados$X - Mu # HQ<-sum(res*drop(solve(SIGMA,res))) ## vamos fazer isto usando a dec. choleski z <- forwardsolve(tchSIGMA,res) HQ <- sum(z*z) detSIGMA <- 2*sum(log(diag(tchSIGMA))) ## det. calculado a partir da choleski if(conc) llik <- -0.5 * (dados$n * (log(2*pi) + log(HQ/dados$n)) + detSIGMA + dados$n) else llik<- -0.5 * (dados$n * log(2*pi) + detSIGMA + HQ) return(llik) } debug(ll) llik <- ll(pars, ...) if(printpars){ print(unlist(pars1)) debug(ll) print(llik) } return(llik) } dados2.var1<-cbind(DADOS2[1:120,1],DADOS2[1:120,2], DADOS2[1:120,3]) dados2.var2<-cbind(DADOS2[c(1:60,121:180),1],DADOS2[c(1:60,121:180),2], DADOS2[c(1:60,121:180),4]) coords2.<-cbind(c(dados2.var1[,1],dados2.var2[,1]), c(dados2.var1[,2],dados2.var2[,2])) var<-c(dados2.var1[,3],dados2.var2[,3]) lista<-list(n=length(var), n1=length(var)/2, n2=length(var)/2, X=var, d=as.vector(dist(coords2.))) ini <-c(a=0.5,rho=0.65,nu1=1,nu2=0.7,tau1=0.4,tau2=sqrt(0.3), sig2=sqrt(0.6),sig1=sqrt(0.65),mu1=1,mu2=3) fixos=NULL debug(geral) geral(ini,pars2=fixos,dados=lista, conc=T) undebug(geral) ajuste <- optim(ini, geral, pars2=fixos, dados=lista,method="BFGS", control=list(fnscale=-1, maxit=3000), conc=T, printpars=F, logp=T) pars.est<-list(a<-ajuste[[1]][[1]], rho<-ajuste[[1]][[2]], nu1<-ajuste[[1]][[3]], nu2<-ajuste[[1]][[4]], tau1<-ajuste[[1]][[5]], tau2<-ajuste[[1]][[6]], sig2<-ajuste[[1]][[7]], sig1<-1 ) names(pars.est) <- c('a','rho','nu1','nu2','tau1','tau2','sig2','sig1') dados<-list(n1=length(var)/2, n2=length(var)/2, X=var, d=as.vector(dist(coords2))) param<-function(pars.est,dados){ SIGMA <- montaSigma(pars.est, dados$d) tchSIGMA <- t(chol(SIGMA)) ONE <- kronecker(diag(1, 2), rep(1, dados$n1)) ONEchSIGMA <- solve(tchSIGMA, ONE) mus <- solve(crossprod(ONEchSIGMA), crossprod(ONEchSIGMA, solve(tchSIGMA, dados$X))) dif<-dados$X-ONE%*%mus sig1<-(1/(dados$n1+dados$n2))*(t(dif)%*%solve(tchSIGMA,dif)) sig2<-pars.est$sig2/sig1 tau1<-pars.est$tau1/sig1 tau2<-pars.est$tau2/sig1 nu1<-pars.est$nu1 nu2<-pars.est$nu2 a<-pars.est$a rho<-pars.est$rho mu1<-mus[1,] mu2<-mus[2,] parametros <- c(mu1,mu2,sig1,sig2,tau1,tau2,nu1,nu2,a,rho) #names(parametros) = c('mu1','mu2','sig1','sig2','tau1','tau2','nu1','nu2','a','rho') return(parametros) } parametros <- param(pars.est,dados)