\documentclass[pnumromarab,normaltoc]{abnt} 

\usepackage[brazil]{babel}
\usepackage[utf8]{inputenc}
\usepackage{abnt-alf}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{alltt}
\usepackage{caption}
\usepackage{textcomp}
\usepackage[T1]{fontenc} 
\usepackage{aurical}    
\usepackage[dvips]{graphicx}
\usepackage{epic,curves}
\usepackage{exscale,relsize}
\usepackage{enumerate}

\makeatletter 
\begin{document}
\noindent \textbf{INFERÊNCIA BAYESIANA}\\
\begin{verbatim}
## carregando pacotes
require(geoR) # necessário para usar rinvchisq
require(geoComp)
require(MCMCpack)
require(MASS)
## Lendo dados
dados <- pivo[,c(6,7,8,1,2)]
## Transformando os dados para a classe geoComp
dados <- as.geoComp(dados)
estima.bayes <- infBayes(dados,tune=c(0.5,0.5,0.5,0.5,0.5),
                         queima=1000,salto=10,nsim=12000,mediahiper=c())
write.table(estima.bayes[[1]],'estima.bayes1.txt')
write.table(estima.bayes[[2]],'estima.bayes2.txt')
## Opção para o usuário entrar com os parâmetros das prioris:
#estima.bayes <- infBayes(dados,tune=c(0.5,0.5,0.5,0.5,0.5),
                          queima=1000,salto=10,nsim=12000,
                          mediahiper=c(1,1,1,1,1,1,1,1,1,1))
summary(estima.bayes[[2]])
res.mh <- resul.mh(resultado=estima.bayes)
res.mh.ic <- data.frame(res.mh[[2]])
names(res.mh.ic) <- c('2.5%','media','97.5%')
write.table(res.mh.ic,'res.mh.ic.txt')
metro.bayes <- data.frame(res.mh[[1]])
write.table(metro.bayes,'metro.bayes.txt')
bor <- cbind(c(0,seq(0,200,l=100),0),
             c(0,sqrt(200^2-seq(0,200,l=100)^2),0))
gr <- pred_grid(bor, by=4)
cokri.bayes <- cokri.bayes(esti.par=metro.bayes,locations=gr,
                           dados.comp=dados)
write.table(cokri.bayes[[1]],'cokri1.txt')
write.table(cokri.bayes[[2]],'cokri2.txt')

## Fazendo o mapa
tes1 <- read.table('cokri1_by4ns1200.txt')#composição
tes2 <- read.table('cokri2_by4ns1200.txt')#coordenadas
tes1 <- as.matrix(tes1)
tes2 <-as.matrix(tes2)
par(mfrow=c(1,3), mar=c(3,3,.5,.5), mgp=c(1.7,0.7, 0),cex=0.9)
xleg <- c(85, 198); yleg <- c(185, 200)
image(structure(list(predict=tes1[,1]), class="kriging"),loc=tes2,bor,
      col=bpy.colors(15),x.leg=xleg, y.leg=yleg,ylim=c(0,220))
text(x=165,y=170,"(A)Areia/IB")
image(structure(list(predict=tes1[,2]), class="kriging"),loc=tes2,bor,
      col=bpy.colors(15),x.leg=xleg, y.leg=yleg,ylim=c(0,220))
text(x=165,y=170,"(B)Silte/IB")
image(structure(list(predict=tes1[,3]), class="kriging"),loc=tes2,bor,
      col=bpy.colors(15),x.leg=xleg, y.leg=yleg,ylim=c(0,220))
text(x=165,y=170,"(C)Argila/IB")

### Fazendo diagnóstico
require(coda)
metro.emv <- read.table('metro.bayes_by4ns1200.txt')
metro.emv <- as.mcmc(metro.emv)
metro.emv.GR <- read.table('metro.bayes.GR.txt')
metro.emv.GR <- as.mcmc(metro.emv.GR)
fac <- autocorr.diag(metro.emv)
crosscor <-crosscorr(metro.emv)
hpd <-HPDinterval(metro.emv)
plot(density(metro.emv[,1]), main='post.mu1')
rug(metro.emv[,1])
lines(seq(unique(hpd[,1][1]),unique(hpd[,2][1]),l=10),rep(0,10),
      type="l",lty=1)
teste <- geweke.diag(metro.emv,frac1=0.90,frac2=0.5)
## Rodando uma segunda cadeia com 
## theta1.inicial <-c(1,1,1,40,0.8) na FunçãoInfBayes.R
## para fazer o teste Gelman eRubin
metro.emv <- read.table('mbmh1.txt')
metro.emv.GR<- read.table('mbmh2.txt')
metro.emv<- as.mcmc(metro.emv)
metro.emv.GR<- as.mcmc(metro.emv.GR)
dadoPgelman <- list()
dadoPgelman[[1]] <- metro.emv
dadoPgelman[[2]] <- metro.emv.GR
res.gelman <- gelman.diag(dadoPgelman,confidence=0.95,transform=FALSE,
                          autoburnin=FALSE)
gelmanplot <- gelman.plot(dadoPgelman,bin.width=10,max.bins=50,
                          confidence=0.95, transform=FALSE,
                          auto.layout=TRUE)
\end{verbatim}
\end{document}
