# Programa_artigo rm(list=ls(all=TRUE)) require(geoR) require(RandomFields) miles2km <- 1.608 data(weather) #157pontos head(weather) (sdP <- sd(weather[, 1])) (sdT <- sd(weather[, 2])) weather[,1]<-weather[, 1] / sdP weather[,2]<-weather[, 2] / sdT head(weather) ind=seq(1:nrow(weather)) set.seed(29) amostra<-sample(ind,52,rep=FALSE) ind2<-ind ind2[amostra]<-0 weather[ind2,] ##################################################################### ####Para predição: 52pontos coords.pred<-weather[amostra, 3:4] var1.pred<-weather[amostra, 1] var2.pred<-weather[amostra, 2] PRED<-cbind(coords.pred,var1.pred,var2.pred) #plot das coordenadas(preto) e retiradas(vermelho) todas<-weather[,3:4] retirados<-weather[amostra, 3:4] plot(todas, xlim=c(-130,-113),ylim=c(39,53)) par(new=TRUE) plot(retirados,col="red",axes=FALSE, ann=FALSE, pch=16, xlim=c(-130,-113),ylim=c(39,53)) ################################ ###1)Dados co-locados: Utilizados na análise coords1<-weather[ind2,3:4] vars1<-weather[ind2,1:2] DADOS1<- cbind(coords1,vars1) 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çã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) } 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 if(logpars){ pars$sig1 <- exp(pars$sig1) pars$sig2 <- exp(pars$sig2) pars$tau1 <- exp(pars$tau1) pars$tau2 <- exp(pars$tau2) #pars$rho <- exp(pars$rho)-1/(1+exp(pars$rho)) RODANDO COM ISSO DA ERRO } if(conc){ pars$sig1 <- 1 } if(with(pars, abs(rho)>1||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))) #Chutes iniciais ini <-c(tau1=sqrt(15),tau2=sqrt(0.1),rho=-0.45,nu1=1.3,nu2=1.5,sig2=sqrt(2),a=50) fixos <- c(mu1=0.5,mu2=0.2,sig1=sqrt(1)) #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) ajuste pars.est<-list(a<-ajuste[[1]][[7]], rho<- (exp(ajuste[[1]][[3]])-1)/(1+exp(ajuste[[1]][[3]])), nu1<-ajuste[[1]][[4]], nu2<-ajuste[[1]][[5]], tau1<-(ajuste[[1]][[1]]), tau2<-(ajuste[[1]][[2]]), sig2<-(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)) sig1<-sig1. #sig2<-(pars.est$sig2/sig1.)^2 sig2<-pars.est$sig2 #tau1<-(pars.est$tau1/sig1.)^2 tau1<-pars.est$tau1 #tau2<-(pars.est$tau2/sig1.)^2 tau2<-pars.est$tau2 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') parametros ################## ###Para Predição:# ################## k = NULL k = as.data.frame(list(mu1,mu2,sig1,sig2,tau1,tau2,nu1,nu2,a,rho)) names(k) = c('mu1','mu2','sig1','sig2','tau1','tau2','nu1','nu2','a','rho') d = as.vector(dist(coords1)) S<-montaSigma(k, d) Sinv<-solve(S) object<-as.list(c(n1=length(vars1[,1]), n2=length(vars1[,2]), mu1=parametros$mu1, mu2=parametros$mu2, sigma1=parametros$sig1, sigma2=parametros$sig2, phi=parametros$a, kappa1=parametros$nu1, kappa2=parametros$nu2 )) locations<-coords.pred geodata1<-as.geodata(DADOS1,coords.col=1:2,data.col=3) geodata2<-as.geodata(DADOS1,coords.col=1:2,data.col=4) pred<-function (object, locations,Sinv, variable.to.predict = 1) { call.fc <- match.call() if (!is.vector(variable.to.predict) || !is.numeric(variable.to.predict) || length(variable.to.predict) != 1) stop("variável deve ser um escalar igual a 1 ou 2") locations <- locations geodata1 <- geodata1 #geodata1 <- get(attr(object, "geodata1")) geodata2 <- geodata2 #geodata2 <- get(attr(object, "geodata2")) #dimnames(locations) <- list(NULL, NULL) n0 <- nrow(locations) n1 <- object$n1 n2 <- object$n2 n <- n1 + n2 res12 <- c((geodata1$data - object$mu1), (geodata2$data - object$mu2)) S012 <- matrix(0, nrow = length(res12), ncol = nrow(locations)) if (any(variable.to.predict == 1)) { mu0 <- rep(object$mu1, n0) var0 <- sum(object$sigma1) S012[1:n1, ] <- cov.spatial(loccoords(geodata1$coords, locations), cov.pars = c((object$sigma1)^2, object$phi), cov.model = "matern", kappa = object$kappa1) S012[(n1 + 1):(n1 + n2), ] <- cov.spatial(loccoords(geodata2$coords, locations), cov.pars = c(rho * sig1 * sig2, object$phi), cov.model = "matern", kappa = (object$kappa1+object$kappa2)/2) } if (any(variable.to.predict == 2)) { mu0 <- rep(object$mu2, n0) var0 <- rep(object$sigma2,n0) S012[1:n1, ] <- cov.spatial(loccoords(geodata1$coords, locations), cov.pars = c(rho * sig1 * sig2, object$phi1), cov.model = "matern", kappa = (object$kappa1+object$kappa2)/2) S012[(n1 + 1):(n1 + n2), ] <- cov.spatial(loccoords(geodata2$coords, locations), cov.pars = c((object$sigma2)^2, object$phi), cov.model = "matern", kappa = object$kappa2) } res <- list() res$predict <- mu0 + drop(crossprod(S012, Sinv) %*% res12) res$krige.var <- var0 - geoR:::.diagquadraticformXAX(S012, lowerA = Sinv[lower.tri(Sinv)], diagA = diag(Sinv)) res$call <- call.fc attr(res, "prediction.locations") <- call.fc$locations return(res) } predição1<-pred(object,locations,Sinv,variable.to.predict=1) predição2<-pred(object,locations,Sinv,variable.to.predict=2) ###Predição! gr <- expand.grid(seq(0,1, by=0.05), seq(0,1, by=0.02)) #grid para predi??o ##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