################################################ ####### Leitura dos dados ######################### setwd("C:\\Users\\Lisa\\Documents\\Arquivos Elisângela\\FMRP- USP2\\Doutorado\\2 Semestre 2012\\Geoestatística\\Artigo") dados<-read.csv2("Dados3.csv",sep=";",dec=",") attach(dados) ###################################### ### Passando para geo dados############## require(geoR) dG <- as.geodata(dados, coords.col=c(1,2), data.col=3) class(dG) ### Ilustração############################ plot(dG) plot(dG, low=T) points(dG, ylab="Coordenadas em Y(UTM)", xlab="Coordenadas em X (UTM)") points(dG, pt.div="quint", cex.min=1, cex.max=1, main="Visualização dos pontos") ###################################### #### Verificando transformação box-cox#### require(MASS) boxcox(PbB~coord.x+ coord.y) # Observa-se que a variável não necessita de transformação ##### Obtendo o variograma ########## v <- variog(dG, max.dist=838544,pairs.min = 15) plot(v, main="Semivariograma: Níveis de Chumbo") ###### Obtendo envelope simulado ############ v.mc <- variog.mc.env(dG, obj.variog=v) plot(v, env=v.mc) ####### Ajuste no olho ############## ef <- eyefit(v) #### Visualizando o grido onde serão feitas as predições#### gr <- expand.grid(seq(826938,838543, len=31),seq(7652011,7661928, len=31)) points(dG) points(gr, pch=19, cex=0.25, col=2) ## Fazendo a predição espacial (krigagem) ########## kc1 <- krige.conv(dG, loc=gr, krige=krige.control(obj.m=ef[[1]]) #Exponencial kc2 <- krige.conv(dG, loc=gr, krige=krige.control(obj.m=ef[[2]]) )#Esférico kc3 <- krige.conv(dG, loc=gr, krige=krige.control(obj.m=ef[[3]]) )#Gaussiano ### Mapa da krigagem########33 image(kc1, col=terrain.colors(21), main=" Mapa de krigagem:Níveis de chumbo") #escala cinza image(kc1, col=gray(seq(1,0, l=21))) ###Mapa de krigagem e contorno image(kc1, col=terrain.colors(21), main=" Mapa de krigagem e contorno:Níveis de chumbo") contour(kc1, add=T, nlev=21) #escala cinza image(kc1, col=gray(seq(1,0, l=21)), main=" Mapa de krigagem e contorno:Níveis de chumbo") contour(kc1, add=T, nlev=21) ## Mapa em perspectiva persp(kc1, theta=20, phi=35, main=" Mapa em perspectiva:Níveis de Chumbo") ######################################################### ######## Abordagem bayesiana ############################## ## Definindo as prioris MC<-model.control() PC<-prior.control(phi.discrete=seq(0,5,len=41),phi.prior="rec")# os outros estão com prioris #default OC<-output.control(n.post=1000) ### Outras prioris MCc <- model.control() PGC <- prior.control(beta.prior="flat", sigmasq.prior="rec", phi.prior="rec", phi=3000, phi.discrete=seq(0,1000,by=20), tausq.rel=0.1) OC <- output.control(n.post=2000) ####Obtendo a krigagem bayesiana ####### set.seed(288901) kb2<-krige.bayes(dG,model=MCc,loc=gr,prior=PGC,output=OC) #### Mapas da krigagem bayesiana ######### image(kb2) image(kb2, col=terrain.colors(21)) image(kb2, col=gray(seq(1,0, l=21))) image(kb2, col=gray(seq(0,1, l=4))) points(dG, add=T) #####Mapa de contorno ####### contour(kb2) ###Mapa da krigagem +contorno image(kb2, col=terrain.colors(21)) contour(kb2, add=T, nlev=21) ###Mapa em perspetiva persp(kb2, theta=20, phi=35) names(kb2$posterior) amostra = (kb2$posterior$sample) beta0 = amostra$beta; sigmasq = amostra$sigmasq; phi = amostra$phi; tausq = (amostra$tausq.rel)*sigmasq #Estudo de convergencia: nburn = 1001:10000 #"Esquecendo" as 100 primeiras iterações amostra.par = kb2$posterior$sample library(coda) parametros = mcmc(cbind(beta0, sigmasq, phi, tausq)[nburn,]) summary(parametros) #Gerando Gráficos range = 1001:10000 #"Esquecendo" as 1000 primeiras iterações #Para os betas pdf('trace1.pdf') par(mfrow = c(4, 1)) ts.plot(beta0[range], type = "l", main = expression(beta[0]), xlab = "Amostra", col="blue") ts.plot(sigmasq[range],type="l", main = expression(sigma^2), xlab = "Amostra", col="77") ts.plot(phi[range], type="l",main = expression(phi),col="hotpink", xlab="Amostra") ts.plot(tausq[range], type = "l", main = expression(tau^2), xlab = "Amostra", col="orange") dev.off() pdf('cor1.pdf') acf(beta0[range], main = expression(beta[0]),col="blue") dev.off() pdf('cor2.pdf') par(mfrow = c(3, 1)) acf(sigmasq[range], main = expression(sigma^2),col="77") acf(phi[range], main = expression(phi),col="hotpink") acf(tausq[range], main = expression(tau^2),col="orange") dev.off() ## visualizando as posterioris (marginais) ## beta|y hist(kb2$posterior$sample[,1], prob=T) density(kb2$posterior$sample[,1]) lines(density(kb2$posterior$sample[,1])) rug(kb2$posterior$sample[,1]) #Mapa da predição pdf('bayes_dens1.pdf') my.colors1 <- colorRampPalette(c("white", " green","yellow","orange","pink","red")) image(s100.kb, xlab="UTM-X", ylab="UTM-Y", val=kb2$predictive$mean, col=my.colors1(100),xlim=c(299570,300100), ylim=c(7504600,7505220),x.leg=c(299800,300100),y.leg=c(7504605,7504625), main="Níveis de Chumbo") dev.off()