# Programa_artigo rm(list=ls(all=TRUE)) require(geoR) ## Parâmetros segundo a parametrização do artigo de Matern Biv Parcimonioso LP <- list(nu1 = 0.9 , nu2 = 0.6, ## parametros da matern a = 1/4 , ## scala da fc de correlação (1/phi) rho = 0.85 , ## parametro da intensidade de correlação entre variaveis sig1 = sqrt(1.9), sig2 = sqrt(1.5),## parametros de scala do termo espacial tau1 = sqrt(0.4), tau2 = sqrt(0.2) ## parametros de escala do termo não espacial ) set.seed(29) gs <- cbind(runif(340, 0, 25), runif(340, 0, 25))#todas coordenadas dim(gs) plot(gs, asp=1, pch=19, cex=0.5) 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 } ## função para montar matrix de convariancia segundo modelo Mat. Biv. Parc. ## (somente para dados colocados!!!) ## assume dados "empilhados" por variavel [Y1, Y2]' 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ção nos parâmetros não satisfeita") 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) } ## simulando dados segundo o modelo Mat. Biv. Parc. dV <- as.vector(dist(gs)) system.time(SIGMA <- montaSigma(LP, distPares=dV)) cholSIGMA <- chol(SIGMA) var<- drop(crossprod(cholSIGMA, rnorm(nrow(cholSIGMA)))) var length(var) ## adicionando as médias de cada variável var <- var + c(rep(2, 340), rep(10, 340)) ## removendo objetos não mais necessarios rm(dV, SIGMA, cholSIGMA) ##Separando as variáveis var1<-var[1:(length(var)/2)] var1 length(var1) var2<-var[((length(var)/2)+1):length(var)] var2 length(var2) c(2, with(LP, c(sig1^2+tau1^2 , sqrt(sig1^2+tau1^2)))) c(mean(var1), var(var1), sd(var1)) c(10, with(LP, c(sig2^2+tau2^2 , sqrt(sig2^2+tau2^2)))) c(mean(var2), var(var2), sd(var2)) plot(variog(coords=gs, data=var1), ty="b"); abline(h=var(var1),lty=3) lines(variog(coords=gs, data=var2), ty="b", col=2); abline(h=var(var2), col=2,lty=3) ##################################################################### ####Para predição: 100pontos coords.pred<-gs[241:340, 1:2] var1.pred<-var1[241:340] var2.pred<-var2[241:340] PRED<-cbind(coords.pred,var1.pred,var2.pred) ##################################################################### ####Separados: var1<-var1[1:240] var2<-var2[1:240] DADOS<-cbind(gs[1:240, 1:2],var1[1:240],var2[1:240]) ################################ ###1)Dados co-locados: Utilizados na análise DADOS1<- DADOS[1:120,1:4] coords1<-DADOS1[,1:2] vars1<-DADOS1[,3:4] 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$tau1 <=0 || pars$tau2 <=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 #se pelo menos um em pars1 conc=F 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) } var<-c(vars1[,1],vars1[,2]) lista<-list(n=length(var), n1=length(var)/2, n2=length(var)/2, X=var, d=as.vector(dist(coords1))) #USANDO CONC=T ini <-c(tau1=sqrt(0.6),tau2=sqrt(0.4),rho=0.7,nu1=1,nu2=1,sig2=sqrt(2),a=30) fixos <- c(mu1=1,mu2=2,sig1=sqrt(2)) #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]][[7]], rho<-ajuste[[1]][[3]], nu1<-ajuste[[1]][[4]], nu2<-ajuste[[1]][[5]], tau1<-exp(ajuste[[1]][[1]]), tau2<-exp(ajuste[[1]][[2]]), sig2<-exp(ajuste[[1]][[6]]), 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(coords1))) 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/sqrt(sig1))^2 #elevado a 2 e tau1 tau 2 tau1<-(pars.est$tau1/sqrt(sig1))^2 tau2<-(pars.est$tau2/sqrt(sig1))^2 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) parametros<-list(mu1<-parametros[[1]], mu2<-parametros[[2]], sig1<-parametros[[3]], sig2<-parametros[[4]], tau1<-parametros[[5]], tau2<-parametros[[6]], nu1<-parametros[[7]], nu2<-parametros[[8]], a<-parametros[[9]], rho<-parametros[[10]] ) names(parametros) = c('mu1','mu2','sig1','sig2','tau1','tau2','nu1','nu2','a','rho') ################## ###Para Predição:# ################## pred1<-krige.control(type.krige = "ok", trend.d = "cte", trend.l = "cte", beta=mu1, cov.model="mat", cov.pars=c(sqrt(sig1),a), kappa=nu1, nugget=tau1 ) pred2<-krige.control(type.krige = "ok", trend.d = "cte", trend.l = "cte", beta=mu2, cov.model="mat", cov.pars=c(sqrt(sig2),a), kappa=nu2, nugget=tau2 ) var1.geodata<-as.geodata(DADOS1,coords.col=1:2,data.col=3) var2.geodata<-as.geodata(DADOS1,coords.col=1:2,data.col=4) ###Predição! gr <- expand.grid(seq(0,1, by=0.05), seq(0,1, by=0.02)) #grid para predição #predição.var1<-ksline(var1.geodata,loc=var1.geodata$coords,data=var1.geodata$data, #cov.pars=c(sig1,a),nugget=tau1,kappa=nu1) #Variável1 predição1<-krige.conv(var1.geodata, locations=gr,krige=pred1) image(predição1, col=terrain.colors(15)) image(predição1, col=gray(seq(.8, 0, leng=15))) #Variável2 predição2<-krige.conv(var2.geodata, locations=gr,krige=pred2) image(predição2, col=terrain.colors(15)) image(predição2, col=gray(seq(.8, 0, leng=15))) ####Predição nas coordenadas utilizadas! predição.var1<-krige.conv(var1.geodata, locations=var1.geodata$coords,krige=pred1) #image(predição.var1, col=terrain.colors(15)) #image(predição.var1, col=gray(seq(.8, 0, leng=15))) predição.var2<-krige.conv(var2.geodata, locations=DADOS1[,1:2],krige=pred2) #image(predição.var2, col=terrain.colors(15)) #image(predição.var2, col=gray(seq(.8, 0, leng=15))) ####Predição dos valores retirados! var1.gd<-as.geodata(PRED,coords.col=1:2,data.col=3) predição.var11<-krige.conv(var1.gd, locations=PRED[,1:2],krige=pred1) str(predição.var11) var2.gd<-as.geodata(PRED,coords.col=1:2,data.col=4) predição.var22<-krige.conv(var2.gd, locations=PRED[,1:2],krige=pred2) str(predição.var22) plot(predição.var11$pred, var1.pred, asp=1); abline(0,1) plot(predição.var22$pred, var2.pred, asp=1); abline(0,1) ##Loglik do modelo loglik_ loglik_ ##Erro Quadrático Médio EQM_v1=(predição.var11$pred - var1.pred)^2 EQM_v1 summary(EQM_v1) EQM_v2=(predição.var22$pred - var2.pred)^2 EQM_v2 summary(EQM_v2) ##Erro Absoluto EA_v1=predição.var11$pred - var1.pred EA_v1 summary(EA_v1) EA_v2=predição.var22$pred - var2.pred EA_v2 summary(EA_v2) #AIC p= 10#número de parametros AIC_ AIC_ ##Intervalo de Predição #A) Das variáveis utilizadas krige.var1<-predição.var1$krige.var krige.var2<-predição.var2$krige.var krige.var<-c(krige.var1,krige.var2) predict1<-predição.var1$predict predict2<-predição.var2$predict predict<-c(predict1,predict2) I.P_inf.<-predict-1.96*sqrt(krige.var) I.P_inf. I.P_sup.<-predict+1.96*sqrt(krige.var) I.P_sup. INTERVALOO<-c(I.P_inf.) variaveis.<-c(var1.geodata$data,var2.geodata$data) I.O..<-(variaveis.>I.P_inf. & variaveis.I.P_inf._ & variav