\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{ESTUDO DOS INTERVALOS DE COBERTURA:}\\
\begin{verbatim}
## Função para simular um conjunto de dados
simu.dados <- function(theta,n.simul){
  mu1 <- theta[1]
  mu2 <- theta[2]
  s1 <- theta[3]
  s2 <- theta[4]
  tau1 <- theta[5]
  tau2 <- theta[6]
  phi <- theta[7]
  rho <- theta[8]
  ## Gerando erros de uma normal bivariada com mu=0, 
  ## variâncias tau1, tau2 e correlação rho 
  v <- matrix(c(tau1^2,tau1*tau2*rho,tau1*tau2*rho,tau2^2),2,2)
  Sc <- kronecker(diag(n.simul/2),v)
  set.seed(3333)
  e <- mvrnorm(n.simul,mu=c(0,0),v)
  ## Gerando um campo aleatório gaussiano com média 0, variância 1 e 
  ## parâmetros de alcance phi 
  set.seed(3333)
  U <- grf(n.simul,grid='irreg',cov.model='exp',cov.pars=c(1,phi),
           mean=0,messages=FALSE)
  S1 <- s1*U$data
  S2 <- s2*U$data
  ## Calculando e organizando os dados
  y1 <-as.geodata(data.frame(U$coords[,1],U$coords[,2],mu1+S1+e[,1]))
  y2 <- as.geodata(data.frame(U$coords[,1],U$coords[,2],mu2+S2+e[,2]))
  seq1 <- seq(1,length(y1$data)*2,by=2)
  seq2 <- seq(2,length(y2$data)*2,by=2)
  y <- c()
  y[seq1] <- y1$data
  y[seq2] <- y2$data
  coords <- data.frame(U$coords[,1],U$coords[,2])
  ## Preparando a saída
  retorna <- list()
  retorna[1][[1]] <- y1$data
  retorna[2][[1]] <- y2$data
  retorna[3][[1]] <- y
  retorna[4][[1]] <- coords
  return(retorna)
}
\end{verbatim}

\begin{verbatim}
simulando <- function(theta,n.simu=1000,n.amostra=300,
                      metodo="L-BFGS-B"){
  MATRIZ <- matrix(nrow=n.simu,ncol=27) 
  dados <- list()
  assign("dados", dados, envir = .GlobalEnv)
  X <- cbind(rep(1:0,length=n.amostra*2),rep(0:1,length=n.amostra*2))
  assign("X", X, envir = .GlobalEnv)
  y <- c()
  assign("y", y, envir = .GlobalEnv)
  for( i in 1:n.simu){
    sim <- simu.dados(theta,n.simul=n.amostra)# n.simul=n. de
                                              # observações simuladas
    y1     <- sim[1][[1]]
    y2     <- sim[2][[1]]
    y      <- sim[3][[1]]
    assign("y",y, envir= .GlobalEnv)
    coords <- sim[4][[1]]
    dados[1][[1]] <- sim[3][[1]]
    dados[2][[1]] <- sim[4][[1]]
    assign("dados", dados, envir = .GlobalEnv)
    estimado <- try(mec(y1,y2,X,coords=coords,metodo=metodo,
                    otimizador="optim"))
    recurso <- list()
    recurso[[1]] <- data.frame(matrix(ncol=5,nrow=8))
    recurso[[2]] <- c(1,1,1)
    if(class(estimado) == "try-error"){
      est <- recurso
    }
    if(class(estimado) != "try-error"){
      est <- estimado
    }
    print(est)
    t1 <- as.matrix(est[[1]][c(2,4,5)])[1,]
    t2 <- as.matrix(est[[1]][c(2,4,5)])[2,]
    t3 <- as.matrix(est[[1]][c(2,4,5)])[3,]
    t4 <- as.matrix(est[[1]][c(2,4,5)])[4,]
    t5 <- as.matrix(est[[1]][c(2,4,5)])[5,]
    t6 <- as.matrix(est[[1]][c(2,4,5)])[6,]
    t7 <- as.matrix(est[[1]][c(2,4,5)])[7,]
    t8 <- as.matrix(est[[1]][c(2,4,5)])[8,]
    linha <- c(matrix(c(t1,t2,t3,t4,t5,t6,t7,t8,est[[2]][1],est[[2]][2],
               est[[2]][3]),nrow=1))
    MATRIZ[i,] <- linha
   }
   MATRIZ <- data.frame(MATRIZ)
   names(MATRIZ) <- c("mu1","Lmin.mu1","Lmax.mu1","mu2","Lmin.mu2",
                      "Lmax.mu2","s1","Lmin.s1","Lmax.s1","s2" ,
                      "Lmin.s2","Lmax.s2","t1","Lmin.t1","Lmax.t1",
                      "t2","Lmin.t2","Lmax.t2" ,"phi","Lmin.phi",
                      "Lmax.phi","rho","Lmin.rho","Lmax.rho","logvero",
                      "logvero.est","conv")              
   return(MATRIZ)
}
\end{verbatim}

\noindent \underline{\textbf{Script para obtenção dos intervalos de cobertura}}  
\begin{verbatim}
## Carregando pacotes
require(MASS)
require(geoR)
require(geoComp)
## Carregar as funções simu.dados e simulando
theta <- c(50,70,2.5,3.5,2.5,3.5,0.6,0.5)
LBFGSB <- simulando(theta,n.simu=1000,n.amostra=300,metodo="L-BFGS-B")
write.table(LBFGSB,"EstudoSimns1000na300.txt")
\end{verbatim}

\noindent \underline{\textbf{Gráficos de densidades a \textit{posteriori}}} 
\begin{verbatim}
dad.sim.na300ns1000 <- read.table('EstudoSimns1000na300.txt')
par(mfrow=c(2,5),mar=c(3.5,3.5,1.5,1.5),mgp=c(1.6,0.6,0))
hist(dad.sim.na300ns1000$mu1,prob=T,main='',xlim=c(44,56),
     ylim=c(0,0.25),xlab=('mu1 estimado'),ylab=('Densidade'))
lines(density(dad.sim.na300ns1000$mu1),col=2)
hist(dad.sim.na300ns1000$mu2,prob=T,main='',xlim=c(60,80),
     ylim=c(0,0.20),xlab=('mu2 estimado'),ylab=('Densidade'))
lines(density(dad.sim.na300ns1000$mu2),col=2)
hist(dad.sim.na300ns1000$s1,prob=T,main='',xlim=c(0,5),
     ylim=c(0,0.9),xlab=('s1 estimado'),ylab=('Densidade'))
lines(density(dad.sim.na300ns1000$s1),col=2)
hist(dad.sim.na300ns1000$s2,prob=T,main='',xlim=c(0,8),
     ylim=c(0,0.6),xlab=('s2 estimado'),ylab=('Densidade'))
lines(density(dad.sim.na300ns1000$s2),col=2)
hist(dad.sim.na300ns1000$t1,prob=T,main='',xlab=('t1 estimado'),
     ylab=('Densidade'))
lines(density(dad.sim.na300ns1000$t1),col=2)
hist(dad.sim.na300ns1000$t2,prob=T,main='',xlim=c(2,5),
     ylim=c(0,2.5),xlab=('t2 estimado'),ylab=('Densidade'))
lines(density(dad.sim.na300ns1000$t2),col=2)
hist(dad.sim.na300ns1000$phi,prob=T,main='',ylim=c(0,3),
     xlab=('phi estimado'),ylab=('Densidade'))
lines(density(dad.sim.na300ns1000$phi),col=2)
hist(dad.sim.na300ns1000$rho,prob=T,main='',xlim=c(0.2,0.7),
     ylim=c(0,10),xlab=('rho estimado'),ylab=('Densidade'))
lines(density(dad.sim.na300ns1000$rho),col=2)
hist(dad.sim.na300ns1000$logvero,prob=T,main='',xlim=c(1450,1700),
     ylim=c(0,0.02),xlab=('logvero'),ylab=('Densidade'))
lines(density(dad.sim.na300ns1000$logvero),col=2)
hist(dad.sim.na300ns1000$logvero.est,prob=T,main='',xlim=c(-1560,-1420),
     ylim=c(0,0.025),xlab=('logvero estimada'),ylab=('Densidade'))
lines(density(dad.sim.na300ns1000$logvero.est),col=2)
\end{verbatim}

\noindent \underline{\textbf{Boxplots das \textit{posterioris} dos parâmetros}} 
\begin{verbatim}
par(mfrow=c(2,4),mar=c(3.5, 3.5, 1.5, 0.5), mgp=c(1.8,.8, 0))
boxplot(dad.sim.na300ns1000$mu1,main='mu1',ylab="Estimativas")
boxplot(dad.sim.na300ns1000$mu2,main='mu2',ylab="Estimativas")
boxplot(dad.sim.na300ns1000$s1,main='s1',ylab="Estimativas")
boxplot(dad.sim.na300ns1000$s2,main='s2',ylab="Estimativas")
boxplot(dad.sim.na300ns1000$t1,main='tau1',ylab="Estimativas")
boxplot(dad.sim.na300ns1000$t2,main='tau2',ylab="Estimativas")
boxplot(dad.sim.na300ns1000$phi,main='phi',ylab="Estimativas")
boxplot(dad.sim.na300ns1000$rho,main='rho',ylab="Estimativas")
\end{verbatim}
\end{document}

