require(RandomFields) require(geoR) ##Lendo o conjunto de dados umidade <- read.table("umidade.txt", head=F) dados<-cbind(umidade[,4:101]) ##selecionando de 4 em 4 seq<-seq(1,98,4) data<-dados[,seq] ############################## ANÁLISE DESCRITIVA ############## #### série temporal t<-seq(1,25,1) plot(t,data[1,],type="l",ylab="Umidade",xlab="Tempo",ylim=c(10,30)) for (i in 2:40) lines(t,data[i,],type="l") ###análise espacial g1<-as.geodata(cbind(umidade[,2],umidade[,3],data[,1])) g2<-as.geodata(cbind(umidade[,2],umidade[,3],data[,2])) g3<-as.geodata(cbind(umidade[,2],umidade[,3],data[,3])) g4<-as.geodata(cbind(umidade[,2],umidade[,3],data[,4])) g5<-as.geodata(cbind(umidade[,2],umidade[,3],data[,5])) g6<-as.geodata(cbind(umidade[,2],umidade[,3],data[,6])) g7<-as.geodata(cbind(umidade[,2],umidade[,3],data[,7])) g8<-as.geodata(cbind(umidade[,2],umidade[,3],data[,8])) g9<-as.geodata(cbind(umidade[,2],umidade[,3],data[,9])) g10<-as.geodata(cbind(umidade[,2],umidade[,3],data[,10])) g11<-as.geodata(cbind(umidade[,2],umidade[,3],data[,11])) g12<-as.geodata(cbind(umidade[,2],umidade[,3],data[,12])) g13<-as.geodata(cbind(umidade[,2],umidade[,3],data[,13])) g14<-as.geodata(cbind(umidade[,2],umidade[,3],data[,14])) g15<-as.geodata(cbind(umidade[,2],umidade[,3],data[,15])) g16<-as.geodata(cbind(umidade[,2],umidade[,3],data[,16])) g17<-as.geodata(cbind(umidade[,2],umidade[,3],data[,17])) g18<-as.geodata(cbind(umidade[,2],umidade[,3],data[,18])) g19<-as.geodata(cbind(umidade[,2],umidade[,3],data[,19])) g20<-as.geodata(cbind(umidade[,2],umidade[,3],data[,20])) g21<-as.geodata(cbind(umidade[,2],umidade[,3],data[,21])) g22<-as.geodata(cbind(umidade[,2],umidade[,3],data[,22])) g23<-as.geodata(cbind(umidade[,2],umidade[,3],data[,23])) g24<-as.geodata(cbind(umidade[,2],umidade[,3],data[,24])) g25<-as.geodata(cbind(umidade[,2],umidade[,3],data[,25])) e1<-likfit(g1,ini=c(0.5,0.5)) e2<-likfit(g2,ini=c(0.5,0.5)) e3<-likfit(g3,ini=c(0.5,0.5)) e4<-likfit(g4,ini=c(0.5,0.5)) e5<-likfit(g5,ini=c(0.5,0.5)) e6<-likfit(g6,ini=c(0.5,0.5)) e7<-likfit(g7,ini=c(0.5,0.5)) e8<-likfit(g8,ini=c(0.5,0.5)) e9<-likfit(g9,ini=c(0.5,0.5)) e10<-likfit(g10,ini=c(0.5,0.5)) e11<-likfit(g11,ini=c(0.5,0.5)) e12<-likfit(g12,ini=c(0.5,0.5)) e13<-likfit(g13,ini=c(0.5,0.5)) e14<-likfit(g14,ini=c(0.5,0.5)) e15<-likfit(g15,ini=c(0.5,0.5)) e16<-likfit(g16,ini=c(0.5,0.5)) e17<-likfit(g17,ini=c(0.5,0.5)) e18<-likfit(g18,ini=c(0.5,0.5)) e19<-likfit(g19,ini=c(0.5,0.5)) e20<-likfit(g20,ini=c(0.5,0.5)) e21<-likfit(g21,ini=c(0.5,0.5)) e22<-likfit(g22,ini=c(0.5,0.5)) e23<-likfit(g23,ini=c(0.5,0.5)) e24<-likfit(g24,ini=c(0.5,0.5)) e25<-likfit(g25,ini=c(0.5,0.5)) plot(c(0,15),c(0,2.6),ty="n",xlab="Distância",ylab="Semivariância") lines(e1) lines(e2) lines(e3) lines(e4) lines(e5) lines(e6) lines(e7) lines(e8) lines(e9) lines(e10) lines(e11) lines(e12) lines(e13) lines(e14) lines(e15) lines(e16) lines(e17) lines(e18) lines(e19) lines(e20) lines(e21) lines(e22) lines(e23) lines(e24) lines(e25) par(mfrow=c(1,4)) tx<-seq(1,25,1) beta<-c(e1$beta,e2$beta,e3$beta,e4$beta,e5$beta,e6$beta,e7$beta,e8$beta,e9$beta,e10$beta, e11$beta,e12$beta,e13$beta,e14$beta,e15$beta,e16$beta,e17$beta,e18$beta,e19$beta,e20$beta, e21$beta,e22$beta,e23$beta,e24$beta,e25$beta) plot(tx,beta,type="b",ylab="Média",xlab="Tempo"); abline(h=mean(beta)) tausq<-c(0,0.0232,0.0023,0,0,0.0033,0.2186,0.3527,0.3492,1.0393, 0.3117,0,0.0542,1.1332,0,0.3367,0,0,0.2142,0.2127,0.2309,0.8556,0,0.0056,0.5551) plot(tx,tausq,type="b",ylab="Efeito Pepita",xlab="Tempo"); abline(h=mean(tausq)) sigmasq<-c(e1$sigmasq,e2$sigmasq,e3$sigmasq,e4$sigmasq,e5$sigmasq,e6$sigmasq,e7$sigmasq,e8$sigmasq,e9$sigmasq,e10$sigmasq, e11$sigmasq,e12$sigmasq,e13$sigmasq,e14$sigmasq,e15$sigmasq,e16$sigmasq,e17$sigmasq,e18$sigmasq,e19$sigmasq,e20$sigmasq, e21$sigmasq,e22$sigmasq,e23$sigmasq,e24$sigmasq,e25$sigmasq) plot(tx,sigmasq,type="b",ylab="Variância",xlab="Tempo"); abline(h=mean(sigmasq)) phi<-c(e1$phi,e2$phi,e3$phi,e4$phi,e5$phi,e6$phi,e7$phi,e8$phi,e9$phi,e10$phi, e11$phi,e12$phi,e13$phi,e14$phi,e15$phi,e16$phi,e17$phi,e18$phi,e19$phi,e20$phi, e21$phi,e22$phi,e23$phi,e24$phi,e25$phi) plot(tx,phi,type="b",ylab="Escala",xlab="Tempo"); abline(h=mean(phi)) ######################### ESTIMAÇÃO ESPAÇO-TEMPORAL ################# ####Preparando o conjunto de dados Hx=umidade[,2] Hy=umidade[,3] TT=c(1,25,1) data=data p26<-dados[,98] TOT<-rbind( cbind(Hx,Hy,data[,1],rep(1,40)), cbind(Hx,Hy,data[,2],rep(2,40)), cbind(Hx,Hy,data[,3],rep(3,40)), cbind(Hx,Hy,data[,4],rep(4,40)), cbind(Hx,Hy,data[,5],rep(5,40)), cbind(Hx,Hy,data[,6],rep(6,40)), cbind(Hx,Hy,data[,7],rep(7,40)), cbind(Hx,Hy,data[,8],rep(8,40)), cbind(Hx,Hy,data[,9],rep(9,40)), cbind(Hx,Hy,data[,10],rep(10,40)), cbind(Hx,Hy,data[,11],rep(11,40)), cbind(Hx,Hy,data[,12],rep(12,40)),cbind(Hx,Hy,data[,13],rep(13,40)), cbind(Hx,Hy,data[,14],rep(14,40)), cbind(Hx,Hy,data[,15],rep(15,40)), cbind(Hx,Hy,data[,16],rep(16,40)), cbind(Hx,Hy,data[,17],rep(17,40)), cbind(Hx,Hy,data[,18],rep(18,40)), cbind(Hx,Hy,data[,19],rep(19,40)), cbind(Hx,Hy,data[,20],rep(20,40)),cbind(Hx,Hy,data[,21],rep(21,40)), cbind(Hx,Hy,data[,22],rep(22,40)), cbind(Hx,Hy,data[,23],rep(23,40)), cbind(Hx,Hy,data[,24],rep(24,40)), cbind(Hx,Hy,data[,25],rep(25,40)) ) ####determinando os valores que não precisam ser estimados. f<- function(param){ param[7] <- param[3] param[25] <- param[3] param[29] <- param[3] param[33] <- param[23] param[35] <- param[13] param[36] <- param[14] param[45] <- param[23] return(param) } ####ESTIMAÇÃO ( Modelo Separável ( a separação está na matriz de anisotropia) ) SEP.M<- list( list( model="stable", var=NA , k= list(NA), aniso=list( NA, 0, 0, 0, NaN,0,0,0,0) ) ,"*", list( model="gencauchy", var=1 ,k= list(NA,NA),aniso= list( 0, 0, 0, 0, 0, 0, 0, 0, NA) ) ,"+", list( model="nugget",var=NA ,aniso=list(NaN, 0, 0, 0, NaN, 0, 0, 0, NaN) ),"*", list( model="gencauchy" , var=1 ,k=list(NaN,NaN) , aniso= list( 0, 0, 0, 0, 0, 0, 0, 0, NaN ) ) ) fitnnS = fitvario( x = cbind(Hx,Hy) , T=TT , data = TOT[,3] , model = SEP.M, cross.me = NULL,table=T,transform=f) fitnnS #### Valores já estimados pelo modelo separável que também serão utilizados no modelo não separável var1<-44.12002 ; k1<-0.84504; a1<-0.001310 k2<-1.6449; k3<-0.1293; a9<-0.756 var2<-0.2246 ######Estimação não-separável... k = c( a= 1 , phi=1 , c = k2 , d=NA , psi=1 , dim=2 ) NSST = list( list( model="nsst",var=var1,k=k , aniso = list(a1, 0, 0, 0, a1, 0, 0, 0, a9) ) ,"*", list( model="gencauchy" , var=1 , k= c(k2,k3), aniso=list( 0, 0, 0, 0, 0, 0, 0, 0, a9) ) ,"+", list( model="nugget",var=var2 ,aniso = list(a1, 0, 0, 0, a1, 0, 0, 0, a9) ),"*", list( model="gencauchy" , var=1 ,k= c(k2,k3), aniso= list( 0, 0, 0, 0, 0, 0, 0, 0, a9) ) ) fitNSST = fitvario( x = cbind(Hx,Hy) , T=TT , data = TOT[,3] , model = NSST, cross.me = NULL,table=T) fitNSST #####predição (fixado os valores do parâmetro d (0, 0.25, 0.5, 0.75, 0.95), os valores preditos serão comparados as volares do tempo 26)....... #####d=0 SEP.M<- list( list( model="stable", var=var1 , k= list(k1), aniso=list( a1, 0, 0, 0, a1,0,0,0,0) ) ,"*", list( model="gencauchy", var=1 ,k= list(k2,k3),aniso= list( 0, 0, 0, 0, 0, 0, 0, 0, a9) ) ,"+", list( model="nugget",var=var2 ,aniso=list(a1, 0, 0, 0, a1, 0, 0, 0, a9) ),"*", list( model="gencauchy" , var=1 ,k=list(k2,k3) , aniso= list( 0, 0, 0, 0, 0, 0, 0, 0, a9 ) ) ) i=1 Pred.SEP=rep(0,40) for(i in 1:40) { kk = as.matrix(expand.grid(Hx[i],Hy[i],c(26))) z <- Kriging( krige.method = "S" , x = kk , grid = FALSE , model=SEP.M , given = TOT[,c(1,2,4)] , data =TOT[,3] ) Pred.SEP[i]= z } Pred.SEP R.SEP<-mean((p26-Pred.SEP)^2) M.SEP<-mean(Pred.SEP) V.SEP<-var(Pred.SEP) R.SEP;M.SEP;V.SEP ######d=0.25 k0.25 = c( a= 1 , phi=1 , c = k2 , d=0.25 , psi=1 , dim=2 ) NSST0.25 = list( list( model="nsst",var=var1,k=k0.25 , aniso = list(a1, 0, 0, 0, a1, 0, 0, 0, a9) ) ,"*", list( model="gencauchy" , var=1 , k= c(k2,k3), aniso=list( 0, 0, 0, 0, 0, 0, 0, 0, a9) ) ,"+", list( model="nugget",var=var2 ,aniso = list(a1, 0, 0, 0, a1, 0, 0, 0, a9) ),"*", list( model="gencauchy" , var=1 ,k= c(k2,k3), aniso= list( 0, 0, 0, 0, 0, 0, 0, 0, a9) ) ) i=1 PredN0.25=rep(0,40) for(i in 1:40) { kk = as.matrix(expand.grid(Hx[i],Hy[i],c(26))) z <- Kriging( krige.method = "S" , x = kk , grid = FALSE , model=NSST0.25 , given = TOT[,c(1,2,4)] , data = TOT[,3] ) # simeia tou sample #simulated values PredN0.25[i]= z } PredN0.25 R.0.25<-mean((p26-PredN0.25)^2) M.0.25<-mean(PredN0.25) V.0.25<-var(PredN0.25) R.0.25;M.0.25;V.0.25 ######d=0.5 k0.5 = c( a= 1 , phi=1 , c = k2 , d=0.5 , psi=1 , dim=2 ) NSST0.5 = list( list( model="nsst",var=var1,k=k0.5 , aniso = list(a1, 0, 0, 0, a1, 0, 0, 0, a9) ) ,"*", list( model="gencauchy" , var=1 , k= c(k2,k3), aniso=list( 0, 0, 0, 0, 0, 0, 0, 0, a9) ) ,"+", list( model="nugget",var=var2 ,aniso = list(a1, 0, 0, 0, a1, 0, 0, 0, a9) ),"*", list( model="gencauchy" , var=1 ,k= c(k2,k3), aniso= list( 0, 0, 0, 0, 0, 0, 0, 0, a9) ) ) i=1 PredN0.5=rep(0,40) for(i in 1:40) { kk = as.matrix(expand.grid(Hx[i],Hy[i],c(26))) z <- Kriging( krige.method = "S" , x = kk , grid = FALSE , model=NSST0.5 , given = TOT[,c(1,2,4)] , data = TOT[,3] ) # simeia tou sample #simulated values PredN0.5[i]= z } PredN0.5 R.0.5<-mean((p26-PredN0.5)^2) M.0.5<-mean(PredN0.5) V.0.5<-var(PredN0.5) R.0.5;M.0.5;V.0.5 ######d=0.75 k0.75 = c( a= 1 , phi=1 , c = k2 , d=0.75 , psi=1 , dim=2 ) NSST0.75 = list( list( model="nsst",var=var1,k=k0.75 , aniso = list(a1, 0, 0, 0, a1, 0, 0, 0, a9) ) ,"*", list( model="gencauchy" , var=1 , k= c(k2,k3), aniso=list( 0, 0, 0, 0, 0, 0, 0, 0, a9) ) ,"+", list( model="nugget",var=var2 ,aniso = list(a1, 0, 0, 0, a1, 0, 0, 0, a9) ),"*", list( model="gencauchy" , var=1 ,k= c(k2,k3), aniso= list( 0, 0, 0, 0, 0, 0, 0, 0, a9) ) ) i=1 PredN0.75=rep(0,40) for(i in 1:40) { kk = as.matrix(expand.grid(Hx[i],Hy[i],c(26))) z <- Kriging( krige.method = "S" , x = kk , grid = FALSE , model=NSST0.75 , given = TOT[,c(1,2,4)] , data = TOT[,3] ) # simeia tou sample #simulated values PredN0.75[i]= z } PredN0.75 R.0.75<-mean((p26-PredN0.75)^2) M.0.75<-mean(PredN0.75) V.0.75<-var(PredN0.75) R.0.75;M.0.75;V.0.75 ######d=0.95 k0.95 = c( a= 1 , phi=1 , c = k2 , d=0.95 , psi=1 , dim=2 ) NSST0.95 = list( list( model="nsst",var=var1,k=k0.95 , aniso = list(a1, 0, 0, 0, a1, 0, 0, 0, a9) ) ,"*", list( model="gencauchy" , var=1 , k= c(k2,k3), aniso=list( 0, 0, 0, 0, 0, 0, 0, 0, a9) ) ,"+", list( model="nugget",var=var2 ,aniso = list(a1, 0, 0, 0, a1, 0, 0, 0, a9) ),"*", list( model="gencauchy" , var=1 ,k= c(k2,k3), aniso= list( 0, 0, 0, 0, 0, 0, 0, 0, a9) ) ) i=1 PredN0.95=rep(0,40) for(i in 1:40) { kk = as.matrix(expand.grid(Hx[i],Hy[i],c(26))) z <- Kriging( krige.method = "S" , x = kk , grid = FALSE , model=NSST0.95 , given = TOT[,c(1,2,4)] , data = TOT[,3] ) # simeia tou sample #simulated values PredN0.95[i]= z } PredN0.95 R.0.95<-mean((p26-PredN0.95)^2) M.0.95<-mean(PredN0.95) V.0.95<-var(PredN0.95) R.0.95;M.0.95;V.0.95 ########Predição total ( esta predição é realizada em todo a região considerada, mas considerando os modelos com valores diferentes do parâmetro d.... p1<-seq(-3,10,13/100) p2<-seq(-2,78,80/100) P1<-c(-3,10,13/100) P2<-c(-2,78,80/100) Tk<-c(26,26,1) ##### Krigagem Pred1 <- Kriging( krige.method = "S" , x = P2, y = P1, T=Tk , grid = TRUE, gridtriple=TRUE , model=SEP.M , given = TOT[,c(1,2,4)] , data = TOT[,3] ) Pred2 <- Kriging( krige.method = "S" , x = P2, y = P1, T=Tk , grid = TRUE, gridtriple=TRUE , model=NSST0.25 , given = TOT[,c(1,2,4)] , data = TOT[,3] ) Pred3 <- Kriging( krige.method = "S" , x = P2, y = P1, T=Tk , grid = TRUE, gridtriple=TRUE , model=NSST0.5 , given = TOT[,c(1,2,4)] , data = TOT[,3] ) Pred4 <- Kriging( krige.method = "S" , x = P2, y = P1, T=Tk , grid = TRUE, gridtriple=TRUE , model=NSST0.75 , given = TOT[,c(1,2,4)] , data = TOT[,3] ) Pred5 <- Kriging( krige.method = "S" , x = P2, y = P1, T=Tk , grid = TRUE, gridtriple=TRUE , model=NSST0.95 , given = TOT[,c(1,2,4)] , data = TOT[,3] ) plot(p26,col="red",type="l",ylab="Predições",xlab="Lacalizações",ylim=c(10,23)) lines(Pred.SEP,col=1) lines(PredN0.25,col=3) lines(PredN0.5,col=4) lines(PredN0.75,col=5) lines(PredN0.95,col=7) zlim<-c("real","d=0","d=0.25",d="0.5",d="0.75",d="0.95") legend(32,23,zlim,fill=c(2,1,3,4,5,7), col=c(2,1,3,4,5,7), lty=0.5, cex=0.5) ###imagens Krigagem ####legenda my.legend <- function(lu.x, lu.y, zlim, col, cex=1) { cn <- length(col) filler <- vector("character", length=(cn-3)/2) legend( lu.x , lu.y , y.i=0.03 , x.i=0.1 , legend=c(format(zlim[2], dig=2) , filler , format(mean(zlim), dig=2) , filler , format(zlim[1], dig=2)) , lty=1, col=rev(col) , cex=cex ) } ###krigagens image(p2,p1,asp=1,Pred1[,,1],col=rainbow(100),zlim=range(Pred1), main="d=0") points(Hx,Hy) col<-rainbow(100) zlim<-range(Pred1[,,1]) my.legend(60,30,zlim,col,0.5) image(p2,p1,asp=1,Pred2[,,1],col=rainbow(100),zlim=range(Pred2), main="d=0.25") image(p2,p1,asp=1,Pred3[,,1],col=rainbow(100),zlim=range(Pred3), main="d=0.5") image(p2,p1,asp=1,Pred4[,,1],col=rainbow(100),zlim=range(Pred4), main="d=0.75") image(p2,p1,asp=1,Pred5[,,1],col=rainbow(100),zlim=range(Pred5), main="d=0.95")