% todas as figuras geradas pelo R sao jogadas no diretorio figs
\documentclass[10pt,twoside]{article}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{RME_p}
\usepackage[brazil,brazilian]{babel}
\usepackage[utf8]{inputenc}

\usepackage{graphicx}   % para inclusão de figuras

 \unespjournalvolume{20}
 \unespjournalnumber{1}
 \unespjournalyear{2000}
 \unespjournalpages{1-10}
 \setcounter{page}{1}

\usepackage{Sweave}
\begin{document}

\title{Um Modelo Geoestatístico Bivariado para Dados Composicionais}
\author{ANA BEATRIZ TOZZO MARTINS\thanks{%
PPGMNE/DES, CESEC/CCE, UFPR/UEM, Rua Mariano Torres, 714, Ap64, CEP: $80060-120$, Curitiba,
Pr, Brazil, E-mail: {\it abtmartins@uem.br}} \and
PAULO JUSTINIANO RIBEIRO JUNIOR\thanks{%
DEst, LEG, UFPR, Caixa Postal: 19.081 , CEP: 81531-990, Curitiba,
Pr, Brazil, E-mail: {\it paulojus@leg.ufpr.br}} \and 
WAGNER HUGO BONAT\thanks{PPGMNE, LEG, UFPR, Caixa Postal: 19.081, CEP: 81531-990, Curitiba, Pr, Brazil, E-mail: {\it wbonat@leg.ufpr.br}} } \maketitle

\begin{resumo}
Este trabalho é motivado pelo interesse em modelar o padrão espacial de dados composicionais. A categoria de problemas de interesse envolve, por exemplo, as frações granulométricas de um solo ou composição química de uma rocha e é composta por dados
 onde as observações são partes de algum todo e são referenciadas
 espacialmente. O interesse central está em predizer conjuntamente o
 padrão espacial dos componentes na área de estudo respeitando a
 restrição da soma dos componentes corresponderem ao total.
Neste sentido, combina-se a teoria de dados composicionais originalmente
 desenvolvida para observações independentes com métodos geoestatísticos
 multivariados. A abordagem proposta declara explicitamente um modelo
 paramétrico que descreve a espacialização das variáveis.
 Em particular neste trabalho, o objetivo é propor e implementar um modelo
 geoestatístico bivariado para dados composicionais obtendo inferências
 via verossimilhança e expressões para predições espaciais, apresentando,
 discutindo e disponibilizando resultados e implementação computacional.
 A metodologia proposta é utilizada na análise de um conjunto de dados
 simulados, cujas composições são formadas por três componentes, possibilitando a obtenção de mapas de predição para
 cada um dos componentes.
\end{resumo}

\begin{palavra:chave}
geoestatística multivariada; dados composicionais; verossimilhança.
\end{palavra:chave}

\section{Introdução}

\par Este estudo é motivado pelo interesse em modelar e descrever o padrão espacial de dados composicionais. Neste sentido, combina-se a teoria de dados composicionais originalmente desenvolvida para observações  independentes (AITCHISON, $1986$) com métodos geoestatísticos (PAWLOWSKY-GLAHN; OLEA, $2004$), adotando-se a declaração explícita do modelo que descreve a espacialização das variáveis (DIGGLE; RIBEIRO JR, $2007$).

\par Dados composicionais consistem de vetores, denominados composições, cujos componentes $X_1,...,X_B$ representam frações de algum ``todo'', e satisfazem a restrição de que a soma dos componentes é igual a um (AITCHISON, $1986$), ou seja,\\
$$X_1 \geq 0,X_2 \geq 0,...,X_B \geq 0,
\hspace{1cm}
\text{e}
\hspace{1cm}
X_1 + X_2 +\cdots+X_B =1. $$
\par O espaço amostral é o simplex unitário de dimensão igual ao número de componentes dado por
$$ \mathbb{S}^B = \{ \b{X} \in \mathbb{R}^B;\,\, X_i > 0, \,\, i=1,...,B;\,\, \b{1}'\b{X}=1 \}.$$
\par Um vetor $\b{W}$ cujos componentes são positivos e medidos na mesma escala denomina-se base. Uma base pode se tornar uma composição através do operador fechamento, $\mathcal{C}$, que garante que a restrição de soma igual a um seja satisfeita:
$$
\begin{array}{cccl}
\mathcal{C} \ : & \! \mathbb{R}^B_+ & \! \longrightarrow & \mathbb{S}^B  \\
                & \! \b{W}          & \! \longrightarrow & \mathcal{C}\big(\b{W}\big)=\dfrac{\b{W}}{\b{1}'\b{W}}.
\end{array}
$$
\par Neste espaço amostral, o simplex, as operações matemáticas de soma e multiplicação definidas no espaço real equivalem às operações pertubação 
$$ \b{X}_1 \oplus \b{X}_2 = (X_{11},X_{12},...,X_{1B}) \oplus (X_{21},X_{22},...,X_{2B}) = \mathcal{C}(X_{11}X_{21},...,X_{1B}X_{2B}),$$ 
e potência 
$$ \alpha \odot (X_{11},X_{12},...,X_{1B}) = \mathcal{C}(X_{11}^{\alpha},X_{12}^{\alpha},...,X_{1B}^{\alpha}),$$
respectivamente, e a média passa a ser a média geométrica $g(\b{X}_1)= \sqrt[B]{\prod_{j=1}^B X_{1j}}$.
\par Uma característica desse tipo de dados é que estes apresentam um efeito de correlação espúria. A restrição de que a soma dos componentes deve ser igual a $1$, implica em correlação negativa entre os componentes fazendo com que as correlações não sejam diretamente interpretáveis (GRAF, $2006$), ou seja, as covariâncias estão sujeitas a controles não estocásticos o que implica, segundo Pawlowsky-Glahn e Olea ($2004$), em singularidade da matriz de covariância de uma composição. Com isto, a aplicação de técnicas estatísticas padrão podem levar a resultados inconsistentes. Para contornar este problema, Aitchison ($1986$) propôs, dentre outras, a transformação razão log-aditiva (alr) que generaliza a transformação logística para um vetor composicional de duas partes e é dada por:
$$
\begin{array}{cccl}
\text{alr} \ : & \! \mathbb{S}^B  & \! \longrightarrow &  \mathbb{R}^{B-1}         \\
               & \! \b{X}         & \! \longrightarrow &  \text{alr} \big(\b{X}\big)=\left(\text{ln}\left(\dfrac{X_1}{X_B}\right),\dots,\text{ln}\left(\dfrac{X_{B-1}}{X_B}\right)
                                                                                     \right)'.
\end{array}
$$
Por outro lado, a transformação inversa denominada transformação logística generalizada aditiva (agl) é dada por
$$
\begin{array}{cccl}
\text{agl} \ : & \! \mathbb{R}^{B-1}   & \! \longrightarrow &   \mathbb{S}^B    \\
               & \! \text{alr}(\b{X})  & \! \longrightarrow &  
\text{agl}(\text{alr} (\b{X}) )=\b{X}=  \left( \text{exp} \left( 
                                        \text{ln}\left( \dfrac{X_1}{X_B}\right)\right), \dots,
                                               \text{exp}(0)
                                        \right)'.
\end{array}
$$
\par A representação gráfica de uma amostra de composições pode ser feita através do diagrama ternário, por exemplo no caso em que $B=3$, um triângulo equilátero cujos vértices representam os três componentes da composição (BUTLER, $2008$). De acordo com Graf (2006), considerando que assintoticamente $\b{Y}=\left( \text{ln} (X_1/X_B),\text{ln} (X_2/X_B),...,\text{ln} (X_{B-1}/X_B) \right)'$  tem distribuição Gaussiana $N_{B-1}(\mathbf{\mu}_{\scriptsize L};\,\Sigma_{\scriptsize L})$ em que $\mathbf{\mu}_{\scriptsize L}$ é o vetor de médias das log-razões e $\Sigma_{\scriptsize L}$ a matriz de covariâncias de ordem $(B-1) \times (B-1)$, o domínio de confiança para $\b{Y}$ é limitado pelo elipsóide
$$B_{1-\alpha}'(\b{Y})=\left\{ \b{Y} \in \mathbb{R}^{B-1} | \left(\text{ln}(\b{Y})-\mu_{\scriptsize L} \right)'\Sigma_{\scriptsize L}^{-1}
                       \left(\text{ln}(\b{Y}) -\mu_{\scriptsize L} \right) \leq \chi^2_{B-1;1-\alpha}\right\}$$
em que $\chi^2_{B-1;1-\alpha}$ é o quantil $(1-\alpha)$ da distribuição qui-quadrado com $B-1$ graus de liberdade e o domínio correspondente para $\b{X}=(X_1,X_2,...,X_B)$ é um subconjunto do simplex $\mathbb{S}^B$
$$B_{1-\alpha}'(\b{X})=\left\{ \b{X} \in \mathbb{S}^B | \left(\text{ln}\frac{\b{X}_{-B}}{X_B} -\mu_{\scriptsize L} \right)'\Sigma^{-1}
\left(\text{ln}\frac{\b{X}_{-B}}{X_B} -\mu_{\scriptsize L} \right) \leq \chi^2_{B-1;1-\alpha}\right\}.$$

\par A teoria de dados composicionais desenvolvida por Aitchison ($1986$) para amostras 
independentes, foi extendida por Pawlowsky-Glahn e Olea ($2004$) 
considerando a dependência espacial segundo métodos geoestatísticos
em uma abordagem que evita declarar completa e explicitamente o modelo
multivariado espacial associado às composições. 
Por outro lado, modelagens multivariadas espaciais são descritas em Diggle e Ribeiro Jr ($2007$), 
Schmidt e Gelfand ($2003$), Banerjee, Carlin e Gelfand ($2004$), 
Schmidt e Sansó ($2006$) garantindo, por construção, 
matrizes de covariância definidas positiva. 
Considera-se aqui o modelo Gaussiano bivariado de componentes comum
proposto por Diggle e Ribeiro Jr ($2007$)
e adotado por Bognola et. al. ($2008$) que utilizam uma variável física como informação secundária.
\par O objetivo deste trabalho é estender o uso do modelo geoestatístico bivariado para estruturas de dados composicionais, derivando e implementando estimação baseada na verossimilhança e obtendo preditores espaciais, na escala original dos dados no simplex, que permitam a construção de mapas de predição dos componentes, ou funcionais destes, na área de estudo. 
\medskip

\section{Metodologia}
\par Para $\b{$X$}\hspace{-0.2cm}=\hspace{-0.2cm}(X_1,...,X_B)'$ uma composição com $B$ componentes e $\b{$Y$}\hspace{-0.2cm}~=~\hspace{-0.2cm}\left(\hspace{-0.1cm}\text{ln}\left(\hspace{-0.1cm}\dfrac{X_1}{X_B}\hspace{-0.1cm}\right),\dots,
\text{ln}\left(\hspace{-0.1cm}\dfrac{X_{B-1}}{X_B}\hspace{-0.1cm}\right)\hspace{-0.1cm}\right)'$ 
um vetor com  $B-1$ elementos, o modelo geoestatístico com componente comum pode ser obtido seguindo a formulação dada em Diggle e Ribeiro Jr (2007). Por simplificação, neste trabalho considera-se composições de apenas $3$ componentes, $X_1=\text{Componente 1}$, $X_2=\text{Componente 2}$ e $X_3=\text{Componente 3},$ que resultam em vetores bivariados.

\par A partir do modelo geoestatístico apresentado em Diggle e Ribeiro Jr ($2007$), 
propõe-se uma adaptação e o modelo pode ser escrito como:
\begin{eqnarray*}\left\{
\begin{array}{ccl}
 Y_{1}(\b{$x$}_i) & = &\mu_1(\b{$x$}_i)+\sigma_{1}U(\b{$x$}_i;\phi)+Z_1(\b{$x$}_i) \\
 Y_{2} (\b{$x$}_{i^{'}}) & = & \mu_2(\b{$x$}_{i^{'}})+\sigma_{2}U(\b{$x$}_{i^{'}};\phi)+Z_2(\b{$x$}_{i^{'}}).    \\
\end{array}
\right.
\label{eq:bivarComU}
\end{eqnarray*}
em que $\b{$x$}_i,\b{$x$}_{i'} \in \mathbb{R}^{2}$, são as localizações amostrais $i,i^{'}=1,...,n_1$, onde $n_1$ é o tamanho da amostra;  $Y_1=\text{ln}(X_1/X_3)$, $Y_2=\text{ln}(X_2/X_3)$  são as variáveis resposta do modelo  de modo que  $\b{$Y$}_{n\times 1}=\left(Y_1(\b{$x$}_1),Y_2(\b{$x$}_1),\dots,Y_1(\b{$x$}_{n_{1}}),Y_2(\b{$x$}_{n_{1}})\right)^{'}$. Neste modelo, assume-se que $U$ é um efeito aleatório de média zero e variância unitária 
com distribuição Gaussiana multivariada com correlações espaciais dadas 
 pela função de correlação exponencial ($\rho_U$), 
caracterizada por um parâmetro $\phi$ que controla o decaimento da correlação como função
da separação espacial entre duas localizações. 
No modelo bivariado geral as unidades de medida são preservadas nas constantes padronizadoras 
$\sigma_{1}$ e $\sigma_{2}$, enquanto que no contexto considerado aqui são admensionais. 
Os efeitos aleatórios $Z_{j} \sim N(0;\tau_j^2)$, $j=1,2$ 
capturam a variabilidade não espacial incluindo a correlação ($\rho$) induzida pela
estrutura composicional.
\par Sendo assim, $\b{$Y$}\sim N_2(\b{$\mu$};\,\mathbf{\Sigma})$,
com  matriz de covariâncias $\mathbf{\Sigma}$ composta pelos elementos
\begin{eqnarray*}
\begin{array}{rclcrcl}
Cov(Y_1(\b{$x$}_i);Y_1(\b{$x$}_i))   & = & \sigma_1^2 +\tau_1^2  &  & 
Cov(Y_1(\b{$x$}_i);Y_1(\b{$x$}_{i'}))& = &  \sigma_1^2 \rho_{U}(\b{$x$}_i;\b{$x$}_{i'})\,   \\
Cov(Y_2(\b{$x$}_i);Y_2(\b{$x$}_i))   & = & \sigma_2^2 +\tau_2^2  &  & 
Cov(Y_2(\b{$x$}_i);Y_2(\b{$x$}_{i'}))& = &  \sigma_2^2 \rho_{U}(\b{$x$}_i;\b{$x$}_{i'})\,
\end{array}
\label{eq:com1}
\end{eqnarray*}
e 
$$Cov(Y_1(\b{$x$}_i);Y_2(\b{$x$}_{i'})) = \sigma_1\sigma_2I_2(i,i')+\tau_1 \tau_2 I_3(i,i')\,$$
onde \\
\begin{minipage}[b]{0.48 \linewidth}
\begin{eqnarray*}I_2(i,i')=\left\{
\begin{array}{lll}
1                                &,\,\,\text{se} & i=i'           \\
\rho_{U}(\b{$x$}_i;\b{$x$}_{i'}) &,\,\,\text{se} & i\neq i'       \\
\end{array}
                           \right.
\label{eq:propostaCom}
\end{eqnarray*}
\end{minipage}\hfill
\begin{minipage}[b]{0.48 \linewidth}
\begin{eqnarray*}I_3(i,i')=\left\{
\begin{array}{lll}
\rho &,\,\,\text{se} & i=i'       \\
   0 &,\,\,\text{se} & i\neq i',  \\
\end{array}
\right.
\label{eq:propostaCom1}
\end{eqnarray*}
\end{minipage}
\noindent Desta forma, a inferência sobre o vetor de parâmetros\\
 $\b{$\theta$}=(\mu_1,\mu_2,\sigma_1,\sigma_2,\tau_1,\tau_2,\phi,\rho)'$ é feita usando a teoria da verossimilhança cuja função de verossimilhança é
$$ L(\b{$\theta$},\b{$Y$})=(2\pi)^{-n/2} |\mathbf{\Sigma}|^{-1/2} \text{exp}
                        \left\{-\frac{1}{2} \left(\b{$Y$}- \b{$\mu$}_{{\scriptsize \b{$Y$}}} \right)'
                        \mathbf{\Sigma}^{-1}
                        \left(\b{$Y$}- \b{$\mu$}_{{\scriptsize \b{$Y$}}} \right)
                                            \right\}.$$
Fazendo-se a reparametrização: $\eta=\sigma_2/\sigma_1$;\,\,\,$\nu_1=\tau_1/\sigma_1$;\,\, $\nu_2=\tau_2/\sigma_1$, pode-se escrever 
$$\mathbf{\Sigma}~=~\sigma_1^2 \textbf{R}~+~\tau_1^2 \textbf{I}_b~=~\sigma_1^2 \textbf{V}$$
em que $\textbf{R}$ corresponde a parte espacial do modelo, $\textbf{I}_b$ é uma matriz bloco diagonal correspondente a parte composicional e
$$
\textbf{V}={\scriptsize \hspace{-0.2cm}
\left[\hspace{-0.2cm}
\begin{array}{ccccccc}
1+\nu_1^2                                   \hspace{-0.2cm}&
\hspace{-0.2cm}\eta+\nu_1\nu_2\rho          \hspace{-0.2cm}&
\rho_{U}(\b{$x$}_1,\b{$x$}_2)               \hspace{-0.2cm}&
\eta\rho_{U}(\b{$x$}_1,\b{$x$}_2)           \hspace{-0.2cm}&
\hspace{-0.1cm} \cdots                      \hspace{-0.1cm}&
\rho_{U}(\b{$x$}_1,\b{$x$}_n)               \hspace{-0.2cm}&
\eta\rho_{U}(\b{$x$}_1,\b{$x$}_n)           \hspace{-0.2cm}  \\

\eta+\nu_1\nu_2\rho                         \hspace{-0.2cm}&
\eta^2+\nu_2^2                              \hspace{-0.2cm}&
\eta\rho_{U}(\b{$x$}_1,\b{$x$}_2)           \hspace{-0.2cm}&
\eta^2 \rho_{U}(\b{$x$}_1,\b{$x$}_2)        \hspace{-0.2cm}&
\hspace{-0.1cm}\cdots                       \hspace{-0.1cm}&
\eta\rho_{U}(\b{$x$}_1,\b{$x$}_n)           \hspace{-0.2cm}& 
\eta^2 \rho_{U}(\b{$x$}_1,\b{$x$}_n)        \hspace{-0.2cm} \\

\rho_{U}(\b{$x$}_2,\b{$x$}_1)               \hspace{-0.2cm}&
\eta \rho_{U}(\b{$x$}_2,\b{$x$}_1)          \hspace{-0.2cm}&
1+\nu_1^2                                   \hspace{-0.2cm}&
\hspace{-0.2cm}\eta+\nu_1\nu_2\rho          \hspace{-0.2cm}&
\hspace{-0.1cm} \cdots                      \hspace{-0.1cm}& 
\rho_{U}(\b{$x$}_2,\b{$x$}_n)               \hspace{-0.2cm}&
\eta \rho_{U}(\b{$x$}_2,\b{$x$}_n)          \hspace{-0.2cm}\\

\eta\rho_{U}(\b{$x$}_2,\b{$x$}_1)           \hspace{-0.2cm}&
\eta^2\rho_{U}(\b{$x$}_2,\b{$x$}_1)         \hspace{-0.2cm}&
\hspace{-0.2cm}\eta+\nu_1\nu_2\rho          \hspace{-0.2cm}& 
\eta^2+\nu_2^2                              \hspace{-0.2cm}&
\hspace{-0.1cm}\cdots                       \hspace{-0.1cm}&
\rho_{U}(\b{$x$}_2,\b{$x$}_n)               \hspace{-0.1cm}&
\eta^2 \rho_{U}(\b{$x$}_2,\b{$x$}_n)        \hspace{-0.2cm} \\

\vdots                                      \hspace{-0.2cm}&
\vdots                                                     &
\vdots                                                     &
\vdots                                      \hspace{-0.2cm}&
\hspace{-0.1cm}\ddots                       \hspace{-0.1cm}&
\hspace{-0.2cm}\vdots                       \hspace{-0.2cm}&
\hspace{-0.2cm}\vdots                       \hspace{-0.2cm} \\

\rho_{U}(\b{$x$}_n,\b{$x$}_1)               \hspace{-0.2cm}&
\eta\rho_{U}(\b{$x$}_n,\b{$x$}_1)           \hspace{-0.2cm}&
\rho_{U}(\b{$x$}_n,\b{$x$}_2)               \hspace{-0.2cm}&
\eta\rho_{U}(\b{$x$}_n,\b{$x$}_2)           \hspace{-0.2cm}&
\hspace{-0.1cm} \cdots                      \hspace{-0.1cm}&
1+\nu_1^2                                   \hspace{-0.2cm}&
\hspace{-0.2cm}\eta+\nu_1\nu_2\rho          \hspace{-0.2cm}\\

\eta \rho_{U}(\b{$x$}_n,\b{$x$}_1)          \hspace{-0.2cm}& 
\eta^2 \rho_{U}(\b{$x$}_n,\b{$x$}_1)        \hspace{-0.2cm}& 
\eta \rho_{U}(\b{$x$}_n,\b{$x$}_2)          \hspace{-0.2cm}& 
\eta^2 \rho_{U}(\b{$x$}_n,\b{$x$}_2)        \hspace{-0.2cm}& 
\hspace{-0.1cm}\cdots                       \hspace{-0.1cm}& 
\eta+\nu_1\nu_2\rho                         \hspace{-0.2cm}& 
\eta^2 \nu_2^2                              \hspace{-0.2cm} \\
\end{array}
\right].
}
$$

\par A função de log-verossimilhança reparametrizada é dada por
% equação (se houver)
\begin{equation}\label{equ:logvero}
l(\b{$\theta$},\b{$Y$})=-\frac{1}{2}\left(n\,\text{ln}(2\pi)+2n\, \text{ln}(\sigma_1) + 
                                             \text{ln}(|\textbf{V}|)+\frac{1}{\sigma_1^2}Qe 
                                    \right).
\end{equation}
\par Considerando $\b{$\mu$}_{{\scriptsize \b{$Y$}}}=\textbf{D}\b{$\mu$}$\hspace{0.2cm} em que \textbf{D} é a matriz do delineamento de ordem $n \times 2$, tem-se 
$$Qe=(\b{$Y$}-\b{$\mu$}_{{\scriptsize \b{$Y$}}})'
     \textbf{V}^{-1}
     (\b{$Y$}-\b{$\mu$}_{{\scriptsize \b{$Y$}}})=\b{$Y$}'\textbf{V}^{-1}\b{$Y$}
                                                 -2(\b{$Y$}'\textbf{V}^{-1}\textbf{D})\b{$\mu$}
                                                +\b{$\mu$}'(\b{$Y$}'\textbf{V}^{-1}\textbf{D})\b{$\mu$}.$$
\par Expressões analíticas fechadas podem ser obtidas para os estimadores de máxima verossimilhança de $\b{$\mu$}=(\mu_1,\mu_2)'$ e $\sigma_1$ diferenciando a função (\ref{equ:logvero}) em relação aos respectivos parâmetros e estes são dados por
\begin{equation}\label{equ:esti}
\begin{array}{rclcrcl}
\hat{\b{$\mu$}} &=& (\textbf{D}'\textbf{V}^{-1}\textbf{D})^{-1}(\textbf{D}'\textbf{V}^{-1}\b{Y}) & \hspace{0.3cm} \text{e} \hspace{0.3cm} & \hat{\sigma}_1 &=& \sqrt{\hat{Q}e/n}.
\end{array}
\end{equation}
É importante notar que $\hat{Q}e$ pode ser escrita como
$$\hat{Q}e=\b{$Y$}'\textbf{V}^{-1}\b{$Y$}
          -(\b{$Y$}'\textbf{V}^{-1}\textbf{D})
           (\textbf{D}'\textbf{V}^{-1}\textbf{D})^{-1}
           (\textbf{D}'\textbf{V}^{-1}\b{$Y$}).$$
Ao substituir as expressões (\ref{equ:esti}) em (\ref{equ:logvero}) obtém-se a função de log-verossimilhança concentrada
$$
l(\b{$\theta$}^*,\b{$Y$})=-\frac{1}{2} \left[ \text{ln}(|\textbf{V}|) + n\left(\text{ln}(2\pi) + \text{ln}(\hat{Qe}) - \text{ln}(n)+ 1\right) \right],
$$
que é uma função do vetor de parâmetros desconhecidos $\b{$\theta$}^*=(\eta,\nu_1,\nu_2,\phi,\rho)'$, e pode ser maximizada numericamente. 
\par Os algoritmos de otimização testados no processo de maximização foram ``L-BFGS-B'', método de Byrd et al. (1995) que permite informar os limites inferior e superior de busca no espaço paramétrico; ``Nelder-Mead'', uma implementação do método de Nelder e Mead (1965); ``Gradiente Conjugado'', baseado no  método de Fletcher e Reeves (1964) e ``BFGS'', um método quasi-Newton. Todos encontram-se implementados no ambiente estatístico R (R Development Core Team, 2008).

\par Do processo de maximização obtém-se $\hat{\b{$\theta$}}^*=(\hat{\eta},\hat{\nu_1},\hat{\nu_2},\hat{\phi},\hat{\rho})'$ e as respectivas variâncias através da matriz Hessiana numérica. A matriz Informação de Fisher observada é definida como o negativo da matriz Hessiana e é dada por
$$\textbf{I}_o(\hat{\b{$\theta$}}^*)=-\dfrac{\partial^2l(\hat{\b{$\theta$}}^*)}{\partial\hat{\b{$\theta$}}^*\partial(\hat{\b{$\theta$}}^*)'}$$
Para se obter $\hat{\mu}_1$, $\hat{\mu}_2$, e $\hat{\sigma}_1$, basta substituir $\hat{\b{$\theta$}^*}$ nas equações em (\ref{equ:esti}).
\par Como o interesse está na obtenção de $\hat{\b{$\theta$}}=(\hat{\mu_1},\hat{\mu_2},\hat{\sigma_1},\hat{\sigma_2},\hat{\tau_1},\hat{\tau_2},\hat{\phi},\hat{\rho})'$ e suas respectivas variâncias, o método delta (DEGROOT, 2002, p.284; COX e HINKLEY, 1974, p.302; AZZALINI, 1996, p.73; PAWITAN, 2001, p.226) é aplicado para obter uma aproximação da distribuição de $\hat{\b{$\theta$}}$. Assintoticamente, a distribuição de $\hat{\b{$\theta$}}$\, será aproximadamente multivariada Gaussiana com vetor de médias $\hat{\b{$\theta$}}=g(\hat{\b{$\theta$}}^*)$\, e variância 
$$\text{Var}(\hat{\b{$\theta$}})=\textbf{I}_{\scriptsize E}(\hat{\b{$\theta$}}) 
                            \geq \nabla g(\hat{\b{$\theta$}}^*)'\hspace{0.1cm}
                                 \textbf{I}_{\scriptsize E}(\hat{\b{$\theta$}}^*)^{-1}  \hspace{0.1cm}
                                 \nabla g(\hat{\b{$\theta$}}^*),$$
em que 
$$
\nabla g(\hat{\b{$\theta$}}^*)= \left(        \dfrac{\partial g(\hat{\b{$\theta$}}^*)}{\partial\eta}, \hspace{0.1cm}
\dfrac{\partial g(\hat{\b{$\theta$}}^*)}{\partial\nu_1},\hspace{0.1cm}
\dfrac{\partial g(\hat{\b{$\theta$}}^*)}{\partial\nu_2},\hspace{0.1cm}
\dfrac{\partial g(\hat{\b{$\theta$}}^*)}{\partial\phi},\hspace{0.1cm}
\dfrac{\partial g(\hat{\b{$\theta$}}^*)}{\partial\rho}
\right)'$$
é a função escore ${\sf U}(\hat{\b{$\theta$}}^*)$. Assim, considerando-se $g(\hat{\b{$\theta$}}^*)$ igual a
$\sigma_2=\eta \sigma_1$, $\tau_1=\nu_1 \sigma_1$, $\tau_2=\nu_2 \sigma_1$, $\phi$ e $\rho$, respectivamente, tem-se
$$\nabla g(\hat{\b{$\theta$}}^*)= \left[\begin{array}{ccccc}
                                  \sigma_1 &       0  &       0 & 0 & 0\\
                                         0 & \sigma_1 &       0 & 0 & 0\\
                                         0 &       0  &\sigma_1 & 0 & 0\\
                                         0 &       0  &       0 & 1 & 0\\
                                         0 &       0  &       0 & 0 & 1\\
                                  \end{array}
                                  \right]$$
e a matriz Informação de Fisher esperada para $\hat{\b{$\theta$}}$ baseada nos dados $\b{Y}$ é substituída pela matriz $\textbf{I}_o(\hat{\b{$\theta$}}^*)$, assintoticamente equivalente, de modo que 
$$\text{Var}(\hat{\b{$\theta$}})\geq \nabla g(\hat{\b{$\theta$}}^*)'\hspace{0.1cm}
                                     \textbf{I}_o(\hat{\b{$\theta$}}^*)^{-1}  \hspace{0.1cm}
                                     \nabla g(\hat{\b{$\theta$}}^*).$$

\par Para encontrar as variâncias para $\hat{\b{$\mu$}}$ e $\hat{\sigma}_1$, através da função (\ref{equ:logvero}), obtém-se
$$\textbf{I}_o(\b{$\mu$})=-\dfrac{\partial^2 l(\b{$\theta$})}{\partial \b{$\mu$}^2} 
                         =\dfrac{1}{\sigma_1^2}(\textbf{D}'\textbf{V}^{-1}\textbf{D})'
  \hspace{0.4cm}\text{e}\hspace{0.4cm}
  \textbf{I}_o(\sigma_1)=-\dfrac{\partial^2 l(\b{$\theta$})}{\partial \sigma_1^2} 
                        =-\dfrac{n}{\sigma_1}+\dfrac{3Qe}{\sigma_1^3},
$$ 
e portanto,
$$
\begin{array}{ccccl}
\text{Var}(\hat{\b{$\mu$}})     &=& \textbf{I}_o(\hat{\b{$\mu$}})^{-1}
                                &=& \hat{\sigma}_1^2 (\textbf{D}'\hat{\textbf{V}}^{-1}\textbf{D})^{-1}\\
\text{Var}(\hat{\sigma}_1)&=& \textbf{I}_o(\hat{\sigma}_1)^{-1}
                                &=& \dfrac{\hat{\sigma}_1^3}{3\hat{Q}e-n\hat{\sigma}_1}.\\
\end{array}
$$
\par O próximo passo é a realização da predição espacial (cokrigagem) de $\b{$Y$}_0$ em localizações não amostradas $\b{$x$}_0=(\b{$x$}_{10},\b{$x$}_{20},...,\b{$x$}_{n_{2}0})$. O vetor de valores esperados correspondentes às variáveis $Y_1$ e $Y_2$ para todas as localizações de predição e a matriz de covariância são dadas pelo seguinte resultado da distribuição Gaussiana multivariada:

\begin{theorem}\label{thm:chave}
Seja $\b{\sf Y}=(\b{$Y$}_0, \b{$Y$})'$  um vetor bivariado com distribuição Gaussiana multivariada conjunta com vetor de médias $\b{$\mu$}=(\mu_{\scriptsize \b{$Y$}_0}, \mu_{\scriptsize \b{$Y$}})'$ e matriz de covariância
$$
\mathbf{\Sigma}=\left[
\begin{array}{cc}
 \mathbf{\Sigma}_{{\scriptsize \b{$Y$}_0}{\scriptsize \b{$Y$}_0}} & 
 \mathbf{\Sigma}_{{\scriptsize \b{$Y$}_0}{\scriptsize \b{$Y$}}} \\
 \mathbf{\Sigma}_{{\scriptsize \b{$Y$}}{\scriptsize \b{$Y$}_0}} & 
 \mathbf{\Sigma}_{{\scriptsize \b{$Y$}}{\scriptsize \b{$Y$}}}
\end{array}
\right]
$$
isto é, $\b{\sf Y} \sim N(\b{$\mu$};\hspace{0.1cm}\mathbf{\Sigma})$. 
Então, a distribuição condicional de $\b{$Y$}_0$ dado $\b{$Y$}$ é também Gaussiana multivariada,
$$ \b{$Y$}_0|\b{$Y$} \sim N(\b{$\mu$}_{\scriptsize {\b{$Y$}_0|\b{$Y$}}};
   \hspace{0.1cm} \mathbf{\Sigma}_{\scriptsize {\b{$Y$}_0|\b{$Y$}}})$$
onde
$$ \b{$\mu$}_{\scriptsize {\b{$Y$}_0|\b{$Y$}}}=
               \b{$\mu$}_{\scriptsize {\b{$Y$}_0}}+
               \mathbf{\Sigma}_{{\scriptsize \b{$Y$}_0}{\scriptsize \b{$Y$}}}
               \mathbf{\Sigma}_{{\scriptsize \b{$Y$}}{\scriptsize \b{$Y$}}}^{-1}
               (\b{$Y$}-\b{$\mu$}_{\scriptsize {\b{$Y$}}})
\hspace{0.5cm}
\text{e}
\hspace{0.5cm} \mathbf{\Sigma}_{\scriptsize {\b{$Y$}_0|\b{$Y$}}}=
               \mathbf{\Sigma}_{{\scriptsize \b{$Y$}_0}{\scriptsize \b{$Y$}_0}}-
               \mathbf{\Sigma}_{{\scriptsize \b{$Y$}_0}{\scriptsize \b{$Y$}}}
               \mathbf{\Sigma}_{{\scriptsize \b{$Y$}}{\scriptsize \b{$Y$}}}^{-1}
               \mathbf{\Sigma}_{{\scriptsize \b{$Y$}}{\scriptsize \b{$Y$}_0}}.
$$
\end{theorem}
\par Sendo desconhecidos os valores de $\b{$\mu$}_{\scriptsize {\b{$Y$}_0}}$, estes são substituídos pelo vetor de médias estimadas obtidas no processo de otimização. A matriz $\mathbf{\Sigma}$, de onde se extrai as matrizes 
$\mathbf{\Sigma}_{{\scriptsize \b{$Y$}_0}{\scriptsize \b{$Y$}_0}}$,
$\mathbf{\Sigma}_{{\scriptsize \b{$Y$}_0}{\scriptsize \b{$Y$}}}$,
$\mathbf{\Sigma}_{{\scriptsize \b{$Y$}}{\scriptsize \b{$Y$}_0}}$ e
$\mathbf{\Sigma}_{\scriptsize {\b{$Y$}\b{$Y$}}}$, 
é obtida substituindo-se os valores estimados para os outros parâmetros na matriz $\textbf{V}$ multiplicada por $\hat{\sigma}_1^2$.
\par Uma vez que a transformação alr foi aplicada aos dados originais e o procedimento de estimação e cokrigagem foi realizada com os dados transformados em $\mathbb{R}^2$, deve-se fazer a transformação de volta do vetor de médias e da matriz de covariância para o espaço amostral original, o simplex $\mathbb{S}^3$, como descrito em Pawlowsky e Olea (2004).
\par O objetivo é calcular para cada localização uma estimativa de 
\begin{equation}\label{equ:intmu}
\b{$\mu$}_{\scriptsize \b{$X$}}=E(\b{$X$})=\int_{\mathbb{S}^B} \b{$X$}f(\b{$X$})d\b{$X$}
\end{equation}
e
\begin{equation}\label{equ:intvar}
\mathbf{\Sigma}_{\scriptsize \b{$X$}}=Cov(\b{$X$},\b{$X$})=\int_{\mathbb{S}^B} (\b{$X$}- \b{$\mu$}_{\scriptsize \b{$X$}})(\b{$X$}- \b{$\mu$}_{\scriptsize \b{$X$}})'f(\b{$X$})d\b{$X$}.
\end{equation}
\par Sabe-se que
$$f(\b{$Y$})=(2\pi)^ {-\frac{B-1}{2}}
             |\mathbf{\Sigma}_{\scriptsize \b{$Y$}}|^ {-\frac{1}{2}}
             \text{exp}\left\{-\dfrac{1}{2} \left(\b{Y}-\b{$\mu$}_{\scriptsize {\b{$Y$}}} 
                                             \right)' 
                                            \mathbf{\Sigma}_{\scriptsize \b{$Y$}}^{-1}
                                            \left(\b{Y}-\b{$\mu$}_{\scriptsize {\b{$Y$}}}
                                            \right)
                       \right\},
$$
segue uma distribuição Gaussiana multivariada e $\b{$Y$}=\text{alr}(\b{$X$})$. Assim, para voltar a escala original, $\b{$X$}=\text{agl}(\b{$Y$})$, é necessário usar o método do jacobiano (JAMES, 2004) cujo jacobiano da transformação é dado por
\begin{equation}\label{equ:jac}
  J_{\scriptsize {\text{alr}(\b{$X$})}}=\left|\dfrac{\partial \b{$Y$}}{\partial \b{$X$}}\right| = \left( \prod_{i=1}^B X_i \right)^{-1}
\end{equation}
de modo que $f(\b{$X$})$ é dada por
$$f(\b{$X$})=(2\pi)^ {-\frac{B-1}{2}}
             |\mathbf{\Sigma}_{\scriptsize \b{$Y$}}|^ {-\frac{1}{2}}
             \left( \prod_{i=1}^B X_i\hspace{-0.1cm} \right)^{-1}\hspace{-0.1cm}
             \text{exp}\left\{-\dfrac{1}{2} \left(\text{alr}(\b{$X$})-\b{$\mu$}_{\scriptsize {\b{$Y$}}} 
                                             \right)' 
                                            \mathbf{\Sigma}_{\scriptsize \b{$Y$}}^{-1}
                                            \left(\text{alr}(\b{$X$})-\b{$\mu$}_{\scriptsize {\b{$Y$}}}
                                            \right)
                       \right\}.
$$
\par Aitchison (1986), Pawlowsky e Olea (2004) afirmam que para se resolver as integrais (\ref{equ:intmu}) e (\ref{equ:intvar}) estas devem ser expressas como
\begin{equation}\label{equ:intmuR}
\b{$\mu$}_{\scriptsize \b{$X$}}=\int_{\mathbb{R}^{B-1}} g_1(\b{${\sf Z}$})f(-\b{${\sf Z}$}'\b{${\sf Z}$})d\b{${\sf Z}$}
\end{equation}
e
\begin{equation}\label{equ:intvarR}
\mathbf{\Sigma}_{\scriptsize \b{$X$}}=\int_{\mathbb{R}^{B-1}} g_2(\b{${\sf Z}$})f(-\b{${\sf Z}$}'\b{${\sf Z}$})d\b{${\sf Z}$}
\end{equation}
em que $\b{${\sf Z}$}$ é a transformação 
\begin{equation}\label{equ:ZZ}
\b{${\sf Z}$}= \dfrac{\text{alr}(\b{$X$})-\b{$\mu$}_{\scriptsize {\b{$Y$}}}}{\sqrt{2}\textbf{\sf R}'}
=\dfrac{1}{\sqrt{2}}(\textbf{\sf R}')^{-1}(\text{alr}(\b{X})-\b{$\mu$}_{\scriptsize {\b{$Y$}}}),
\end{equation}
com \textbf{\sf R} sendo a decomposição Cholesky de $\mathbf{\Sigma}_{\scriptsize \b{$Y$}}$, matriz triangular superior. As integrais (\ref{equ:intmuR}) e (\ref{equ:intvarR}) podem, então, ser aproximadas pela integração de Gauss-Hermite multivariada de ordem ${\sf k}$:
\begin{equation}\label{equ:igh}
\begin{array}{ccc}
\displaystyle \int_{\mathbb{R}^{B-1}}\hspace{-0.1cm} g(\b{${\sf Z}$})f(-\b{${\sf Z}$}'\b{${\sf Z}$})d\b{${\sf Z}$} \hspace{-0.2cm}&\approx&\hspace{-0.2cm}
\displaystyle\sum_{i_1=1}^{{\sf k}}\displaystyle\sum_{i_2=1}^{{\sf k}}\cdots\displaystyle\sum_{i_{B-1}=1}^{{\sf k}}
\omega_{i_1}\omega_{i_2}\cdots\omega_{i_{B-1}}
g({\sf Z}_{i_1},{\sf Z}_{i_2},...,{\sf Z}_{i_{B-1}}).
\end{array}
\end{equation}
Nesta aproximação, os pesos $\omega_{i_1}\omega_{i_2}\cdots\omega_{i_{B-1}}$ e as abscissas ${\sf Z}_{i_1},{\sf Z}_{i_2},...,{\sf Z}_{i_{B-1}}$ são conhecidas e seus valores podem ser encontrados, por exemplo, em Abramowitz e Stegun (1972, p. 924). Segundo Gammerman (1997) ordens de quadratura de $6$ a $8$ são suficientes. 
\par Explicitando o procedimento, observa-se que a integral (\ref{equ:intmu}) pode ser reescrita como
\begin{eqnarray*}
\b{$\mu$}_{\scriptsize \b{$X$}}\hspace{-0.2cm}=\hspace{-0.2cm}
\int_{\mathbb{S}^B} g(\b{$X$})
\left((2\pi)^ {-\frac{B-1}{2}}
             |\mathbf{\Sigma}_{\scriptsize \b{$Y$}}|^ {-\frac{1}{2}}
             \text{exp}\left\{-\dfrac{1}{2} \left(\text{alr}(\b{X})-\b{$\mu$}_{\scriptsize {\b{$Y$}}} 
                                            \right)' 
                                            \mathbf{\Sigma}_{\scriptsize \b{$Y$}}^{-1}
                                            \left(\text{alr}(\b{X})-\b{$\mu$}_{\scriptsize {\b{$Y$}}}
                                            \right)\right\}
\right)d{\b{$X$}}
\label{eq:intmu1}
\end{eqnarray*}
onde
\begin{equation}\label{equ:gg}
 g(\b{$X$})=\b{$X$}\left( \prod_{i=1}^B X_i \right)^{-1}.
\end{equation}
\par Considere que
$\mathbf{\Sigma}_{\scriptsize \b{$Y$}}=\textbf{\sf R}'\textbf{\sf R}$
e
$\mathbf{\Sigma}_{\scriptsize \b{$Y$}}^{-1}=(\textbf{\sf R}'\textbf{\sf R})^{-1}
                                           =\textbf{\sf R}^{-1}(\textbf{\sf R}')^{-1}
                                           =\textbf{\sf R}^{-1}(\textbf{\sf R}^{-1})'.$
Da transformação (\ref{equ:ZZ}), tem-se
$$\b{$X$}=\text{agl}(\b{$\mu$}_{\scriptsize {\b{$Y$}}}+\sqrt{2}\textbf{\sf R}'\b{${\sf Z}$})= \text{agl}(\b{$Y$})$$
e de (\ref{equ:jac}), 
\begin{equation}\label{equ:par} 
\partial \b{$X$}=\left( \prod_{i=1}^B X_i \right)\partial \b{$Y$}.
\end{equation}
Além disso, 
$$J=\left| \dfrac{\partial \b{$Y$}}{\partial \b{${\sf Z}$}} \right|
              =\left| \dfrac{\partial(\b{$\mu$}_{\scriptsize {\b{$Y$}}}+\sqrt{2}
                              \textbf{\sf R}'\b{${\sf Z}$})}{\partial \b{${\sf Z}$}}\right|
              = \left|\sqrt{2}\textbf{\sf R}'\right|=(\sqrt{2})^{B-1}|\textbf{\sf R}'|\,;$$
sendo $\textbf{\sf R}$ a decomposição cholesky de $\mathbf{\Sigma}_{\scriptsize \b{$Y$}}$,
$$|\textbf{\sf R}'|=|\textbf{\sf R}|=\left|\mathbf{\Sigma}_{\scriptsize \b{$Y$}} \right|^{\frac{1}{2}}\,, $$
e
\begin{equation}\label{equ:parX}
  J=\left| \dfrac{\partial \b{$Y$}}{\partial \b{${\sf Z}$}}\right|=
     2^{\frac{B-1}{2}}\left|\mathbf{\Sigma}_{\scriptsize \b{$Y$}}\right|^{\frac{1}{2}}
\hspace{0.5cm} \Rightarrow \hspace{0.5cm}
\partial \b{$Y$}=2^{\frac{B-1}{2}}\left|\mathbf{\Sigma}_{\scriptsize \b{$Y$}} \right|^{\frac{1}{2}}\partial \b{${\sf Z}$}.
\end{equation}
Substituindo a equação (\ref{equ:parX}) em (\ref{equ:par}) tem-se
\begin{eqnarray*}
  \partial \b{$X$}= \left( \prod_{i=1}^B X_i \right) 2^{\frac{B-1}{2}} \left|\mathbf{\Sigma}_{\scriptsize \b{$Y$}}\right|^{\frac{1}{2}} \partial \b{${\sf Z}$} 
\hspace{0.5cm} \Rightarrow \hspace{0.5cm}
\partial \b{${\sf Z}$}=2^{-\frac{B-1}{2}}\left( \prod_{i=1}^B X_i \right)^{-1}\left|\mathbf{\Sigma}_{\scriptsize \b{$Y$}}\right|^{-\frac{1}{2}}\partial \b{$X$}.
\label{equ:parZ}
\end{eqnarray*}
Por outro lado, de (\ref{equ:gg}) tem-se 
$g\left(\text{agl}(\b{$\mu$}_{\scriptsize {\b{$Y$}}}\hspace{-0.1cm}+\hspace{-0.1cm}\sqrt{2}\textbf{\sf R}'
\b{${\sf Z}$})\hspace{-0.1cm}\right)= \text{agl}\left(\b{$\mu$}_{\scriptsize {\b{$Y$}}}\hspace{-0.1cm}+\hspace{-0.1cm}\sqrt{2}\textbf{\sf R}'\b{${\sf Z}$}  \right)\hspace{-0.1cm}\left( \prod_{i=1}^B X_i \right)^{-1}$ que substituída em (\ref{equ:intmuR}) produz
\begin{equation}\label{equ:intmuuuu}
\b{$\mu$}_{\scriptsize \b{$X$}}=
\int_{-\infty}^{+\infty} \int_{-\infty}^{+\infty}
%\pi^{-\frac{b}{2}} \text{agl}\left(\b{$\mu$}_{\scriptsize {\b{$Y$}}}+
%\sqrt{2}\textbf{\sf R}'\b{${\sf Z}$}\right)
\hspace{0.1cm}g_1(\b{${\sf Z}$})\text{exp}\{\b{${\sf Z}$}'\b{${\sf Z}$} \}
d\b{${\sf Z}$} 
\end{equation}
com $g_1(\b{${\sf Z}$})=\pi^{-\frac{B-1}{2}} \text{agl}\left(\b{$\mu$}_{\scriptsize {\b{$Y$}}}+
\sqrt{2}\textbf{\sf R}'\b{${\sf Z}$}\right).$
Logo, a aproximação de Gauss-Hermite (\ref{equ:igh}) de ordem $3$, por exemplo, para $\b{$\mu$}_{\scriptsize \b{$X$}}$ é
\begin{eqnarray*}
 \begin{array}{ccl}
  \b{$\mu$}_{\scriptsize \b{$X$}} &\approx&
   \sum_{i_1=1}^{\sf 3}\sum_{i_2=1}^{\sf 3}\omega_{i_1} \omega_{i_2}
   g_1({\sf Z}_{i_1},{\sf Z}_{i_2}) 
 \end{array}
\end{eqnarray*}
Lembrando que $\b{$\mu$}_{\scriptsize \b{$X$}}$ representa a média da composição em cada localização, esta é um vetor trivariado. Então na função $g_1$, $\b{$\mu$}_{\scriptsize {\b{$Y$}}}$ é um vetor bivariado extraído de $\b{$\mu$}_{\scriptsize {\b{$Y$}_0|\b{$Y$}}}$, cujo primeiro elemento corresponde à média da cokrigagem para a variável $Y_1$, o segundo para a variável $Y_2$ e $\textbf{\sf R}'$ de ordem $2\times 2$ será a correspondente matriz de variância/covariâncias de $Y_1$, $Y_2$, extraída do bloco diagonal da matriz $\mathbf{\Sigma}_{\scriptsize {\b{$Y$}_0|\b{$Y$}}}$. Assim, em cada localização a resolução de (\ref{equ:intmuuuu}) implica na resolução de $3$ integrais.
\par O mesmo procedimento é aplicado para a obtenção de $\mathbf{\Sigma}_{\scriptsize \b{$X$}}$ dada pela integral (\ref{equ:intvarR}), agora reescrita como
\begin{equation}\label{equ:intvarrrr}
  \b{$\Sigma$}_{\scriptsize \b{$X$}} =
  \int_{-\infty}^{+\infty} \int_{-\infty}^{+\infty}
  g_2(\b{${\sf Z}$})\text{exp}\{\b{${\sf Z}$}'\b{${\sf Z}$} \} d\b{${\sf Z}$} 
\end{equation}
com  
$g_2(\b{${\sf Z}$})=\pi^{-\frac{B-1}{2}}
\left( \text{agl}( 
                  \b{$\mu$}_{\scriptsize {\b{$Y$}}}+\sqrt{2}\textbf{\sf R}'\b{${\sf Z}$}
                  )-\b{$\mu$}_{\scriptsize {\b{$X$}}}
\right)
\left( \text{agl}( 
                  \b{$\mu$}_{\scriptsize {\b{$Y$}}}+\sqrt{2}\textbf{\sf R}'\b{${\sf Z}$}
                  )-\b{$\mu$}_{\scriptsize {\b{$X$}}}
\right)'
$
e, neste caso, sendo $\mathbf{\Sigma}_{\scriptsize {\b{$X$}}}$ uma matriz de ordem $3\times3$, a resolução da integral (\ref{equ:intvarrrr}) para cada localização implica na resolução de $9$ integrais. Portanto, a estimação do vetor de médias e variâncias em $n_2$ localizações implica na resolução de $n_2 \times 12$ integrais.
\medskip
%****************************************
%   RESULTADOS       
%****************************************
\section{Análise de Dados Composicionais Simulados}
\par Como exemplo de aplicação da metodologia proposta analisou-se três conjuntos de dados simulados, cada um com $100$ valores de porcentagens de três componentes, $X_1$, $X_2$ e $X_3$ e as localizações em um quadrado unitário encontram-se representadas na Figura~\ref{fig:localpivo}.
\par O primeiro conjunto, denominado configuração~$1$, foi gerado a partir do vetor de parâmetros $\b{$\theta$}\hspace{-0.1cm}=\hspace{-0.1cm}(\mu_1,\mu_2,\sigma_1,\sigma_2,\tau_1,\tau_2,\phi,\rho)=(-0,2;-0,5;1;1,5;0,3;0,3;0,6;0,9)$. Buscou-se, colocando valores das médias próximos e de mesmo sinal, fazer com que a nuvem de pontos se situasse na parte central do triângulo. O fato de estarem concentrados deu-se pelo alto valor de $\rho$: quanto menor, mais espalhados. O fato das variâncias dos efeitos espaciais serem próximos com $\sigma_1<\sigma_2$ fez com que os pontos se aproximassem mais do vértice $X_2$ do que do $X_1$ e, por último, valores iguais para as variâncias $\tau_1$ e $\tau_2$ justificou os pontos no interior do diagrama.\\
%****************************************
%   localização
%****************************************
\setkeys{Gin}{width=0.5\textwidth} %tamanho da figura em relação ao latex
\begin{figure}[htp!]
\begin{center}
\includegraphics{./figs/-001}
\vspace{-1cm}
\end{center}
\caption{Distribuição das localizações da primeira configuração no quadrado unitário.} \label{fig:localpivo}
\end{figure}
%****************************************
\par Em cada uma destas localizações tem-se os percentuais de $X_1$, $X_2$ e $X_3$ e as $100$ composições representadas por círculos em um diagrama ternário podem ser vistas na Figura~\ref{fig:histDigC1}. Composições mais próximas a um vértice têm altas proporções do componente correspondente àquele vértice, de modo que esta figura indica que as maiores porcentagens na amostra ocorreram para componente $X_3$. Também, pode-se observar, por exemplo, que a variabilidade da razão $X_3$/$X_1$ foi maior que a variabilidade da razão $X_2$/$X_1$. Por outro lado, a variabilidade da razão $X_1$/$X_3$, $X_2$/$X_3$ foram similares. Além disso, uma região de confiança de $4$-desvios-padrão contemplou todas as amostras e o centro da distribuição dado pelo fechamento da média geométrica dos três componentes está representado pelo ponto vermelho no diagrama. Pelos histogramas nota-se que as porcentagens dos três componentes não se aproximam de uma distribuição Gaussiana. 
%****************************************
%    elipses 1
%****************************************
\setkeys{Gin}{width=0.67\textwidth} % ERA 0.75  tamanho da figura em relação ao latex  
\begin{figure}[htp!]
\begin{center}
\includegraphics{./figs/-002}
\vspace{-0.8cm}
\end{center}
\caption{Distribuição de $X_1$, $X_2$ e $X_3$ e diagrama ternário das composições para a primeira configuração.} \label{fig:histDigC1}
\end{figure}
%****************************************
\par Porém, ao fazer a transformação alr nos dados simulados obteve-se as variáveis $Y_1=\text{ln}(X_1/X_3)$ e $Y_2=\text{ln}(X_2/X_3)$ e os dados transformados passaram a seguir distribuição Gaussiana como pode ser visto pela Figura~\ref{fig:histCorr1}. O diagrama de dispersão, por sua vez, mostra pontos alinhados em forma linear crescente indicando uma possível correlação linear positiva. 
%****************************************
%    diagrama dispersao 1
%****************************************
\setkeys{Gin}{width=0.90\textwidth} %tamanho da figura em relação ao latex
\begin{figure}[htp!]
\begin{center}
\includegraphics{./figs/-003}
\vspace{-0.5cm}
\end{center}
\caption{Distribuição das log-razão e correspondente diagrama de dispersão para a primeira configuração.} \label{fig:histCorr1}
\end{figure}
%****************************************
\par Em seguida, fez-se o ajuste do modelo. Dentre os métodos de otimização testados adotou-se o método ``L-BFGS-B'' por não apresentar problemas de convergência e as estimativas obtidas para os parâmetros bem como os respectivos intervalos de confiança calculados pelo método delta são apresentados na Tabela~\ref{tab:est1}. Observa-se, em geral, estimativas não muito próximas aos valores dos parâmetros mas com mesmos sinais e que apenas três de oito intervalos contêm os verdadeiros valores.  
%****************************************
% TABELA 1
%****************************************
\begin{table}
 \caption{Estimativas, erros padrão e intervalos de confiança para a primeira configuração, pelo método delta via método de otimização ``L-BFGS-B''.}
\vspace{-0.5cm}
 \label{tab:est1}
 \[ \begin{tabular}{crrrrr} \hline
   Parâmetro &  Valor       & Estimativa & Erro Padrão &   LI. Delta &  LS. Delta\\ \hline
     $\mu_1$ & -0,2         & -0,9925955 &  1,60663318 &  -2,5992287 &  0,6140377\\
     $\mu_2$ & -0,5         & -1,6890643 &  1,98310148 &  -3,6721658 &  0,2940372\\
  $\sigma_1$ &  1           &  1,1530662 &  0,10190384 &   1,0511623 &  1,2549700\\
  $\sigma_2$ &  1,5         &  1,7581358 &  0,05412534 &   1,7040104 &  1,8122611\\
    $\tau_1$ &  0,3         &  0,4004140 &  0,02972426 &   0,3706897 &  0,4301382\\ 
    $\tau_2$ &  0,3         &  0,4274941 &  0,04045833 &   0,3870357 &  0,4679524\\
      $\phi$ &  0,6         &  0,9530621 &  0,51284603 &   0,4402160 &  1,4659081\\
      $\rho$ &  0,9         &  0,9573298 &  0,02626424 &   0,9310656 &  0,9835941\\ \hline
\end{tabular}\]
\end{table}
%****************************************
\par No processo de otimização os valores iniciais para o vetor de médias foram considerados como as médias dos valores observados para $Y_1$ e $Y_2$. Metade da variância calculada para $Y_1$ foi atribuída para o efeito espacial e a outra metade para o efeito composicional. Da mesma forma, procedeu-se com $Y_2$. O valor incial para $\rho$ foi calculado como o coeficiente de correlação de Pearson e $\phi=\text{min}+0,2(\text{max}-\text{min})$, onde ``min'' e ``max'' foram, respectivamente, a menor e maior distância entre duas localizações. Então, estes valores indicaram uma log-verossimilhança concentrada de $-301,35132$. Após o processo de otimização, a log-verossimilhança maximizada foi $59,64009$.
\par Como resultado da cokrigagem obteve-se, para cada localização, o vetor de médias e a matriz de covariância no espaço $\mathbb{R}^2$.  A volta dos valores preditos para o simplex $\mathbb{S}^3$ foi feita em uma grade constituída de $1156$ pontos usando aproximação de Gauss-Hermite com ordem de quadratura $k=7$. Na Figura~\ref{fig:mapagh1ei} tem-se os mapas dos três componentes conforme a configuração 1 e, assim como revelou o diagrama ternário, as maiores proporções ocorreram para o componente $X_3$. A Figura~\ref{fig:obsXpred1}, por sua vez, mostra boa concordância entre os valores observados e preditos para os três componentes. 
%************************************************************************************************
%% MAPA  1  QUADRATURA GAUSSIANA - SEM ESCALA  k=7 ns=1000
%************************************************************************************************
\setkeys{Gin}{width=0.90\textwidth} % 0.99 tamanho da figura em relação ao latex
\begin{figure}[htp!]
\begin{center}
\includegraphics{./figs/-004}
\vspace{-0.5cm}
\end{center}
\caption{Mapas das porcentagens de $X1$, $X2$ e $X3$ para dados da configuração~$1$.} \label{fig:mapagh1ei}
\end{figure}
%****************************************
%  1  observados x preditos  
%****************************************
\setkeys{Gin}{width=0.90\textwidth} %tamanho da figura em relação ao latex
\begin{figure}[htp!]
\begin{center}
\includegraphics{./figs/-005}
\vspace{-0.5cm}
\end{center}
\caption{Valores observados versus preditos de areia, silte e argila para a configuração $1$.} \label{fig:obsXpred1}
\end{figure}
%****************************************
\par O segundo conjunto de dados, configuração $2$, foi gerado considerando-se o vetor de parâmetros $(\mu_1,\mu_2,\sigma_1,\sigma_2,\tau_1,\tau_2,\phi,\rho)=(1;1;1,2;1,5;0,9;1;0,6;0,5)$. Assim como na configuração 1, em relação às médias, buscou-se fazer com que a nuvem de pontos continuasse na parte central do diagrama; a proximidade dos pontos ao vértice $X_3$ justificada pelo baixo valor de $\phi$ e o espalhamento pelo baixo valor de $\rho$. Também aqui, $\sigma_1<\sigma_2$ fez com que os pontos se aproximassem mais de $X_2$ do que de $X_1$ e a diferença entre $\tau_1$ e $\tau_2$ forçou os pontos a se aproximarem mais do lado que liga os vértices $X_1$ e $X_2$, com maior aproximação para $X_2$ por $\tau_1$ ser o menor valor.
%****************************************
\par Para este vetor de parâmetros a configuração das composições no diagrama ternário da Figura~\ref{fig:histDigC2} apresenta-se mais espalhada em relação à anterior com maiores porcentagens relativas aos componentes $X_2$ e $X_3$. Ao mesmo tempo, observa-se muitas composições com porcentagens muito baixas de $X_3$ o que se confirma pelo respectivo histograma.   
%****************************************
%    elipses 2
%****************************************
\setkeys{Gin}{width=0.67\textwidth} %tamanho da figura em relação ao latex
\begin{figure}[htp!]
\begin{center}
\includegraphics{./figs/-006}
\vspace{-0.8cm} 
\end{center}
\caption{Distribuição de $X_1$, $X_2$ e $X_3$ e diagrama ternário das composições para a segunda configuração.} \label{fig:histDigC2}
\end{figure}
%****************************************
\par A Figura~\ref{fig:histCorr2} mostra que os dados transformados também seguem uma distribuição Gaussiana e a nuvem de pontos no diagrama de dispersão continua de forma linear crescente mas com um espalhamento maior em relação ao da configuração $1$. \\
%****************************************
%    diagrama dispersao 2
%****************************************
\setkeys{Gin}{width=0.90\textwidth} %tamanho da figura em relação ao latex
\begin{figure}[htp!]
\begin{center}
\includegraphics{./figs/-007}
\vspace{-0.5cm}
\end{center}
\caption{Distribuição das log-razão e correspondente diagrama de dispersão para a segunda configuração.} \label{fig:histCorr2}
\end{figure}
%****************************************
\par Na Tabela~\ref{tab:est2} tem-se os resultados do ajuste do modelo para a configuração $2$ e o que se observa é que as estimativas permanecem com mesmos sinais; para as médias não são boas mas os respectivos intervalos de confiança contêm os verdadeiros valores. Somente o intervalo para $\sigma_2$ não conteve o valor do parâmetro apesar de estar próximo ao limite inferior. Os valores das log-verossimilhanças concentrada e otimizada foram iguais a $-330,3923$ e $298,63124$, respectivamente.
%****************************************
% TABELA 2
%****************************************
\begin{table}
 \caption{Estimativas, erros padrão e intervalos de confiança para a segunda configuração, pelo método delta via método de otimização ``L-BFGS-B''.}
\vspace{-0.5cm}
 \label{tab:est2}
 \[ \begin{tabular}{crrrrr} \hline
   Parâmetro &  Valor       & Estimativa & Erro Padrão &   LI. Delta &  LS. Delta\\ \hline
     $\mu_1$ &  1           &  0,3883751 &  1,4633944  & -1,07501928 &  1,8517695\\
     $\mu_2$ &  1           &  0,2844116 &  1,7558488  & -1,47143723 &  2,0402604\\
  $\sigma_1$ &  1,2         &  1,2584047 &  0,1046900  &  1,15371472 &  1,3630947\\
  $\sigma_2$ &  1,5         &  1,8313026 &  0,2648620  &  1,56644065 &  2,0961646\\
    $\tau_1$ &  0,9         &  0,9824042 &  0,1806521  &  0,80175207 &  1,1630563\\
    $\tau_2$ &  1           &  1,0052126 &  0,2132289  &  0,79198365 &  1,2184415\\
      $\phi$ &  0,6         &  0,5236357 &  0,4430230  &  0,08061266 &  0,9666587\\
      $\rho$ &  0,5         &  0,5401549 &  0,1562457  &  0,38390917 &  0,6964006\\ \hline
\end{tabular}\]
\end{table}
%****************************************
\par Na Figura~\ref{fig:mapagh2ei} tem-se os mapas onde se observam poucas mudanças em relação aos mapas obtidos para a configuração $1$ e a Figura~\ref{fig:obsXpred2} mostra que concordância já não é tão boa para $X_1$ quanto para $X_2$ e $X_3$. 
%************************************************************************************************
%% MAPA  2   POR QUADRATURA GAUSSIANA - SEM ESCALA  k=7 ns=1000
%************************************************************************************************
\setkeys{Gin}{width=0.90\textwidth} % 0.99 tamanho da figura em relação ao latex
\begin{figure}[htp!]
\begin{center}
\includegraphics{./figs/-008}
\vspace{-0.5cm}
\end{center}
\caption{Mapas das porcentagens de $X1$, $X2$ e $X3$ para dados da configuração~$2$.} \label{fig:mapagh2ei}
\end{figure}
%****************************************
%  2  observados x preditos  
%****************************************
\setkeys{Gin}{width=0.90\textwidth} %tamanho da figura em relação ao latex
\begin{figure}[htp!]
\begin{center}
\includegraphics{./figs/-009}
\vspace{-0.5cm}
\end{center}
\caption{Valores observados versus preditos de areia, silte e argila para a configuração 2.} \label{fig:obsXpred2}
\end{figure}
%****************************************
\par A última configuração deste estudo caracterizou-se pelo conjunto de parâmetros $(\mu_1,\mu_2,\sigma_1,\sigma_2,\tau_1,\tau_2,\phi,\rho)=(-0,2;-1;0.45;0,13;0,3;0,3;0,6;0,95)$. O baixo valor absoluto de $\mu_1$ fez com que os dados se aproximassem do lado esquerdo do diagrama com proximidade a $X_1$; caso contrário, a aproximação seria para o lado direito. Também neste caso, o baixo valor de $\phi$ aproximou os pontos ao vértice $X_3$ enquanto o alto valor de $\rho$ tornou-os concentrados novamente. Como nesta configuração $\sigma_1>\sigma_2$, os pontos se afastaram do vértice $X_2$ e continuaram no interior do diagrama pelos valores iguais e baixos para $\tau_1$ e $\tau_2$.
%****************************************
\par Como pode ser visto na Figura~\ref{fig:histDigC3} tem-se baixos percentuais do componente $X_2$, e as distribuições de $X_1$, $X_2$ e $X_3$ se aproximam de uma distribuição Gaussiana. 
%****************************************
%    elipses 3
%****************************************
\setkeys{Gin}{width=0.65\textwidth} %tamanho da figura em relação ao latex
\begin{figure}[htp!]
\begin{center}
\includegraphics{./figs/-010}
\vspace{-0.8cm}
\end{center}
\caption{Distribuição de $X_1$, $X_2$ e $X_3$ e diagrama ternário das composições para a terceira configuração.} \label{fig:histDigC3}
\end{figure}
%****************************************
\par De acordo com a Figura~\ref{fig:histCorr3}, os dados transformados continuam apresentando distribuição Gaussiana e os pontos no diagrama de dispersão também apresentam-se em forma linear crescente mas com espalhamento maior ainda em relação aos anteriores.
%****************************************
%    diagrama dispersao 3
%****************************************
\setkeys{Gin}{width=0.90\textwidth} %tamanho da figura em relação ao latex
\begin{figure}[htp!]
\begin{center}
\includegraphics{./figs/-011}
\vspace{-0.5cm}
\end{center}
\caption{Distribuição das log-razão e correspondente diagrama de dispersão para a terceira configuração.} \label{fig:histCorr3}
\end{figure}
%****************************************
\par Pela Tabela~\ref{tab:est3}, também neste caso as estimativas conservaram os mesmos sinais dos valores verdadeiros dos parâmetros e observa-se que apenas os intervalos de confiança para as variâncias dos efeitos aleatórios $Z_1$ e $Z_2$ não contêm o verdadeiro valor. Neste caso, o valor da log-verossimilhança concentrada foi $-92,32910$ e o da log-verossimilhança otimizada foi $14,34055$.
%****************************************
% TABELA 3
%****************************************
\begin{table}
 \caption{Estimativas, erros padrão e intervalos de confiança para a terceira configuração, pelo método delta via método de otimização ``L-BFGS-B''.}
\vspace{-0.5cm}
 \label{tab:est3}
 \[ \begin{tabular}{crrrrr} \hline
   Parâmetro &  Valor       & Estimativa & Erro Padrão &   LI. Delta &  LS. Delta\\ \hline
     $\mu_1$ &  -0,2        & -0,5896071 &  0,69119506 &  -1,2808021 &  0,1015880\\
     $\mu_2$ &  -1          & -1,1146150 &  0,40214445 &  -1,5167595 & -0,7124706\\
  $\sigma_1$ &   0,45       &  0,5305950 &  0,09559040 &   0,4350046 &  0,6261854\\
  $\sigma_2$ &   0,13       &  0,1765225 &  0,05364894 &   0,1228735 &  0,2301714\\
    $\tau_1$ &   0,3        &  0,3488848 &  0,03599187 &   0,3128929 &  0,3848766\\
    $\tau_2$ &   0,3        &  0,3446375 &  0,03074784 &   0,3138897 &  0,3753853\\
      $\phi$ &   0,6        &  0,7403627 &  0,52040659 &   0,2199561 &  1,2607693\\
      $\rho$ &   0,95       &  0,9383066 &  0,03744746 &   0,9008592 &  0,9757541\\ \hline
\end{tabular}\]
\end{table}
%****************************************
\par Na Figura~\ref{fig:mapagh3ei} tem-se os mapas referentes à configuração $3$ podendo-se verificar uma grande mudança das predições de areia e silte em relação às outras configurações. Neste caso, a pior concordância entre os valores observados e os preditos ocorreu para o segundo componente (Figura~\ref{fig:obsXpred3}). 
%************************************************************************************************
%% MAPA  3  QUADRATURA GAUSSIANA - SEM ESCALA  k=7 ns=1000
%************************************************************************************************
\setkeys{Gin}{width=0.90\textwidth} % 0.99 tamanho da figura em relação ao latex
\begin{figure}[htp!]
\begin{center}
\includegraphics{./figs/-012}
\vspace{-0.5cm}
\end{center}
\caption{Mapas das porcentagens de $X1$, $X2$ e $X3$ para dados da configuração~$3$.} \label{fig:mapagh3ei}
\end{figure}
%****************************************
%  3  observados x preditos  
%****************************************
\setkeys{Gin}{width=0.90\textwidth} %tamanho da figura em relação ao latex
\begin{figure}[htp!]
\begin{center}
\includegraphics{./figs/-013}
\vspace{-0.5cm}
\end{center}
\caption{Valores observados versus preditos de areia, silte e argila para a configuração 3.} \label{fig:obsXpred3}
\end{figure}
%****************************************
\par Todo o trabalho foi realizado utilizando recursos de \textit{software} livre em ambiente operacional GNU/Linux; no ambiente estatístico R (R Development Core Team,~$2008$), utilizando o pacote \textit{geoR} (RIBEIRO JR; DIGGLE, $2001$), \textit{compositions} (BOOGAART; TOLOSANA; BREN, $2008$), \textit{statmod} (SMYTH; HU; DUNN, $2009$) e rotinas desenvolvidas especificamente para o desenvolvimento deste trabalho.\\
%\citeonline{mccullach+nelder:1989} escreveram um livro que se tornou a referência básica no tema \cite{mccullach+nelder:1989}.
%O artigo de \citeonline{barreto+demetrio:1995} pode ser citado \textit{inline}
%como anteriormente ou como uma referência %\cite{barreto+demetrio:1995}.
%\citeonline{arnold:1992} menciona é uma citação de um capítulo de livro \cite{arnold:1992}.
%Citações duplas também são possíveis % \cite{bussab:1976,stringer+cullis:2002}
%Conclusões (se houverem)

{\noindent \bf \large Discussão}\vspace{.5cm}
\par Este trabalho possibilitou o estudo do comportamento do modelo em três configurações distintas de dados composicionais. De forma geral, o componente $X_1$ apresentou-se o mais similar dentre as configurações com percentuais de médios para baixos. Com relação a $X_2$, a configuração $2$ evidencia a maior variação nas porcentagens e as maiores porcentagens para $X_3$ apareceram na configuração $1$.

\par Pode-se observar pelos diagramas de dispersão que a configuração $1$ que mais se aproximou de uma forma linear foi a que apresentou as piores estimativas, enquanto que a configuração $3$ que esperaria-se pior não resultou a pior em termos de qualidade de ajuste. As melhores estimativas foram obtidas com os dados da configuração $2$. Os intervalos construídos via aproximação quadrática pelo método delta não foram bons principalmente para os parâmetros de variância e métodos alternativos devem ser avalidados como intervalos baseados em verossimilhança perfilhada e intervalos de credibildade sob o enfoque da inferência bayesiana. 

\par Verificou-se em estudos e análises não reportados aqui que os resultados obtidos considerando-se a transformação de volta do vetor de médias e da matriz de covariância para o simplex $\mathbb{S}^3$ por quadratura de Gauss-Hermite de ordem igual a $7$ não diferiram dos obtidos com ordem $20$. Uma outra maneira de se fazer esta transformação seria através de simulação. Neste caso, em primeiro lugar, em cada localização gerariam-se dados bivariados de uma distribuição Gaussiana padrão obtendo-se uma distribuição de percentuais para cada componente; em seguida faria-se uma decomposição da matriz de covariância obtida por cokrigagem por algum método numérico, por exemplo, fatoração cholesky e aplicaria-se uma transformação linear, produto dos dados bivariados pela matriz cholesky para a obtenção do vetor bivariado resposta. Daí, faria-se a transformação agl neste vetor para resgatar a composição. Este procedimento seria repetido para os $1156$ pontos de predição. É possível que se verifique uma melhora nos resultados à medida que o número de simulações aumente. Desta forma, a distribuição dos valores simulados em cada localização seria obtida possibilitando o cálculo de outros funcionais além da média e variância. 

\par Em relação aos mapas de valores preditos construídos para as três configurações, ressalta-se que estes revelam os resultados já observados no diagrama ternário mas agora sob outra forma. De fato, as duas primeiras se mostram mais semelhantes por componente dado as nuvens de pontos localizadas na parte central do triângulo. As mudanças observadas nos mapas correspondentes à configuração $3$ se justificam pelo deslocamento da nuvem de pontos para o lado esquerdo implicando em baixos valores de $X_2$. Por último, as maiores concordâncias entre valores observados e preditos resultaram da configuração $1$. \\

{\noindent \bf \large Conclusões}\vspace{.5cm}
\par Os procedimentos adotados permitiram a construção de mapas de composições com três componentes por uma metodologia que implicitamente garante a restrição de que as frações somem $1$, não só nos pontos simulados como nos pontos preditos. 
\par O modelo proposto captura variações espaciais, induzidas pelas composições e não estruturadas.
\par A declaração explícita do modelo permite que sejam feitas inferências sobre parâmetros de forma usual. Abre-se ainda a possibilidade de tratamento Bayesiano a fim de se considerar nas predições a incerteza associada à estimação dos parâmetros do modelo. 
\par As análises podem ser expandidas para estimação de outros funcionais que não necessariamente resultem em mapas de componentes médios. 
\par Há a necessidade de se investigar alternativas para computação mais eficiente e considerar outras formas de especificação do modelo multivariado para o caso de maiores números de componentes. \vspace{.5cm}\\

%Agradecimentos (se houverem)
{\noindent \bf \large Agradecimentos}\vspace{.5cm}
\par Agradecimento à CAPES pelo apoio financeiro. Esse trabalho foi parcialmente financiado pela FINEP projeto CT-INFRA/UFPR.%Digite os agradecimentos \vspace{.5cm}

%citação do artigo
\info{MARTINS, A. B. T.; RIBEIRO JR, P. J.; BONAT, W. H. Título em inglês. \textit{Rev.
Mat. Estat.,} São Paulo, v.xx, n.x, p.xx-xx, 2000.}
\medskip

\begin{abstract}
A multivariate model for spatial compositions is presented
and results for inference on model parameters and
spatial prediction of the compositions are derived
including  computational implementation.
The methodology is used for the analysis of soil fractions
obtaining maps of sand, silt and clay contents along an unit square.
\end{abstract}
\begin{keyword}
multivariate geostatistic; compositional data; likelihood.
\end{keyword}

\begin{thebibliography}{99}
\bibitem{ABRAMOWITZ}
{ABRAMOWITZ, M.; STEGUN, I. A. {\it Handbook of mathematical functions with formulas, graphs, and mathematical tables}. Washington: Milton Abramowitz and Irene A. Stegun, 1972. 1046p.}

\bibitem{AITCHISON}
{AITCHISON, J. \textit{The statistical analysis of compositional data}. New Jersey: The Blackburn Press, 1986. 416p.}

\bibitem{AZZALINI}
{AZZALINI, A. \textit{Statistical inference based in the likelihood}. London: Chapman and Hall, 1996. 341p.}

\bibitem{BANERJEE}
{BANERJEE, S.; CARLIN, B. P.; GELFAND, G. E. \textit{ Hierarchical modelling and analysis for spatial data}. Boca Raton: Chapman and Hall, 2004. 452p.}

\bibitem{BOGNOLA}
{BOGNOLA, I. A.; RIBEIRO JR, P. J.; SILVA, E. A. A; LINGNAU, C.; HIGA, A. R.
\textit{Modelagem uni e bivariada da variabilidade espacial de rendimento de pinus taeda l.}.
Floresta, v.38, n.2, p.373$-$385, 2008.}

%% olhar o que é isto!!!
\bibitem{BOOGAART}
{BOOGAART, G. v. d.; TOLOSANA, R.; BREN, M. \textit{compositions: compositional data analysis}. Disponível em:$<$http://www.stat.boogaart.de/compositions$>$. Acesso em 25 maio 2009.}

\bibitem{BUTLER}
{BUTLER, A.; GLASBEY, C. A. \textit{Latent Gaussian model for compositional data with zeros}. Journal of the Royal Statistical Society, Series C, v.57, n.5, p.505$-$520, 2008.}

\bibitem{BYRD}
{BYRD, R. H.; LU, P.; NOCEDAL, J.; ZHU, C. \textit{A limited memory algorithm for bound constrained optimization}. SIAM J. Scientific Computing, v.16, p.1190$-$1208, 1995.}

\bibitem{COX}
{COX, D. R.;HINKLEY, D. V. \textit{Theoretical Statistics}. London: Chapman and Hall, 1974. 511p.}

\bibitem{DEGROOT}
{DEGROOT, M. H.; SCHERVISH, M. J. \textit{Probability and Statistics}. 3.ed. USA:
Addison-Wesley, 2002. 816p.}

\bibitem{DIGGLE}
{DIGGLE, P. J.; RIBEIRO JR, P. J. \textit{Model-based geostatistics}. USA:
Springer Series in Statistics, 2007. 228p.}

\bibitem{FLETCHER}
{FLETCHER, R.; REEVES, C. M. \textit{Function minimization by conjugate gradients}. Computer Journal, v.7, p.148$-$154, 1964.}

\bibitem{GAMMERMAN}
{GAMMERMAN, D. \textit{Markov chain monte carlo}. London: Chapman and Hall, 1997. 245p.}

\bibitem{GRAF}
{GRAF, M. \textit{Precision of compositional data in a stratified two-stage cluster sample: comparison of the swiss earnings structure survey 2002 and 2000}. Survey Research Methods Section, ASA, Session 415: Sample Survey Quality V, p.3066$-$3072, 2006.} 

\bibitem{JAMES}
{JAMES, B. R. \textit{Probabilidade: um curso em nível intermediário}. 3.ed. Rio de Janeiro: IMPA, 2004. 304p.}

\bibitem{NELDER}
{NELDER, J. A., MEAD, R. \textit{A simplex algorithm for function minimization}. Computer Journal, v.$7$, p.308$-$313, $1965$.}

\bibitem{PAWITAN}
{PAWITAN, Y. \textit{In all likelihood: statistical modelling and inference using likelihood}. New York: Oxford University Press, 2001. 528p.}

\bibitem{PAWLOWSKY}
{PAWLOWSKY-GLAHN, V.; OLEA, R. A. \textit{Geostatistical analysis of compositional data}. New York: Oxford University Press, Inc., 2004. 181p.}

\bibitem{R}
{R Development Core Team.\textit{R: a language and environment for statistical computing}. Vienna, Austria, 2009. Disponível em:<http://www.R-project.org>. Acesso em: 25 maio 2009.}

\bibitem{RIBEIRO}
{RIBEIRO JR, P. J.; DIGGLE, P. J. \textit{geoR: a package for geostatistical analysis}. R-NEWS, v.1, n.2, p.15$-$18, 2001. Disponível em:<http://CRAN.R-project.org/doc/Rnews/>. Acesso em 13 junho 2009.}

\bibitem{SCHMIDT}
{SCHMIDT, A. M.; GELFAND, A. E. \textit{A bayesian coregionalization approach for multivariate pollutant data}. Journal of Geophysical Research, v.108, p.10$-$1$-$18$-$8, 2003.}

\bibitem{SCHMIDT}
{SCHMIDT, A. M.; SANSÓ, B. Modelagem bayesiana da estrutura de covariância de processos espaciais e espaço temporais. In: 17 SINAPE e ABE-Associação Brasileira de Estatística, 2006, Caxambu. \textit{Minicurso}, Caxambu: Associação Brasileira de Estatística, 2006.}

\bibitem{SMYTH}
{SMYTH, G.; HU, Y.; DUNN, P. \textit{statmod: statistical modeling}. Disponível em:<http://www.statsci.org/r>. Acesso em 25 maio 2009.}
\end{thebibliography}

%%% a ser usado pelo editor

\unespreceived{01.01.2005}

\unespapproved{01.01.2005}
\end{document}
%%%%%%%%%%%%%% 
%    FIM
%%%%%%%%%%%%%%  

