%\RequirePackage{lineno}
\documentclass[a4paper,12pt]{article}
\usepackage{fancyhdr}
\fancyhf{}
\usepackage[latin1]{inputenc} % Acentuação
\usepackage[brazil]{babel}    % Escrevendo em português
\usepackage{graphicx}
\usepackage{latexsym,amssymb,amsmath,amsfonts} %fontes e simbolos de AMS
\usepackage[left=2.5cm, right=2.5cm, top=2.5cm,bottom=2.5cm]{geometry}
\usepackage{url}
\usepackage{indentfirst}
\usepackage{enumerate}
%\usepackage{comment}
\newcommand{\mat}[1]{\mbox{\boldmath{$#1$}}} 

\begin{document}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\pagestyle{fancy}
\setcounter{page}{1}
\lhead{$64^{\textrm{a}}$ RBras e $18^{\textrm{0}}$ SEAGRO }
\rhead{29 de julho à 02 de agosto de 2019, Cuiabá - MT}
\lfoot{}
\cfoot{\thepage}
\rfoot{}
%\setpagewiselinenumbers

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{center}
{\large {\bf Método de Newton e Gauss-Newton na estimação dos parâmetros de modelo de regressão não linear}}
\vspace{0.5cm}
\end{center}

\begin{center}
 {\bf Edilson Marcelino Silva} \footnote[1]{Programa de Pós-Graduação em Estatística e Experimentação Agropecuária, Departamento de Estatística (DES), Universidade Federal de Lavras (UFLA). e-mail: \it 	edilsonmg3@hotmail.com},  {\bf Ariana Campos Fr\"{u}hauf} \footnote[2]{Programa de Pós-Graduação em Estatística e Experimentação Agropecuária, Departamento de Estatística (DES), Universidade Federal de Lavras (UFLA). e-mail: \it 	arianafruhauf@gmail.com}, {\bf Felipe Augusto Fernandes} \footnote[3]{Programa de Pós-Graduação em Estatística e Experimentação Agropecuária, Departamento de Estatística (DES), Universidade Federal de Lavras (UFLA). e-mail: \it 	fernandesfelipest@gmail.com},
{\bf Gustavo Sérgio de Paula} \footnote[4]{Iniciação científica, Departamento de Estatística (DES), Universidade Federal de Lavras (UFLA). e-mail: \it 	gustavo028@hotmail.com}, {\bf Joel Augusto Muniz}  \footnote[5]{Departamento de Estatística (DES), Universidade Federal de Lavras (UFLA). e-mail: \it 	joamuniz@ufla.br }, {\bf Tales Jesus Fernandes}  \footnote[6]{Departamento de Estatística (DES), Universidade Federal de Lavras (UFLA). e-mail: \it 	tales.jfernandes@ufla.br }
 %{\bf Tales Jesus Fernandes}  \footnote[5]{Departamento de Estatística (DES), Universidade Federal de Lavras (UFLA). e-mail: \it 	tales.jfernandes@ufla.br } \\
\end{center} 

\section{Introdução}

Em estatística, os modelos de regressão consistem em ajustar uma função a um conjunto de $n$ observações. A regressão estuda a relação entre duas ou mais variáveis e incorpora os erros experimentais sendo expressa por:
$${\bf y}=f({\bf X},\mat{\beta})+\mat{\varepsilon}$$
em que, ${\bf y}$ é o vetor com os valores observados (variável dependente); ${\bf X}$ matriz com as covariáveis (variável independente); $f$ é a função que associa as variáveis; $\mat{\beta}$ é o vetor de parâmetros e $\mat{\varepsilon}$ é o vetor de erros experimentais.

Os modelos de regressão são classificados basicamente em modelos lineares e não lineares em relação aos seus parâmetros. Considera-se o modelo de regressão como linear quando as derivadas parciais em relação a cada parâmetro do modelo não dependem de nenhum parâmetro do modelo. Por outro lado, se pelo menos uma derivada parcial depender de algum parâmetro este é classificado como não linear. Não há forma fechada para a solução do sistema de equações normais de um modelo não linear, sendo necessário utilizar métodos numéricos iterativos (DRAPER e SMITH, 2014).

Para a solução aproximada de sistema de equações foi criado por Newton (1643-1727) o denominado ``Método de Newton", um processo que é aplicado iterativamente para obter aproximação da solução, tão precisa quanto se queira. O método de Newton aproxima $f(\mat{\beta})$ em $\mat{\beta}^n$ por sua aproximação em série de Talylor até o termo quadrático, assim:
$$f(\mat{\beta}) \approx f(\mat{\beta^n})+{\bf J}(\mat{\beta^n})'(\mat{\beta} - \mat{\beta}^n)+\dfrac{1}{2}(\mat{\beta} - \mat{\beta}^n)'{\bf H}(\mat{\beta^n})(\mat{\beta} - \mat{\beta}^n)$$
em que ${\bf J}$ é a matriz de derivadas parciais de $f(\mat{\beta})$ em relação a $\mat{\beta}$ denominada matriz Jacobiana e ${\bf H}$ é matriz de segundas derivadas denominada de matriz Hessiana. 

Derivando a função quadrática em relação a $\mat{\beta}$ igualando a zero é obtida a condição para ter o mínimo:

$${\bf J}^n+{\bf H}^n(\mat{\beta} - \mat{\beta}^n)=0$$
ou
 $$\mat{\beta} = \mat{\beta^n} - ({\bf H}^n)^{-1}{\bf J}^n.$$ 
 
 Assim a solução da equação $f(\mat{\beta})$ é aproximada pelo método iterativo de Newton como apresentado abaixo:
$$\mat{\beta}^{n+1} = \mat{\beta^n} - ({\bf H}^n)^{-1}{\bf J}^n$$
(BARD, 1974).

Entre as descobertas mais importantes de Gauss (1777-1855) pode-se destacar o método dos mínimos quadrados. O diretor de um observatório na Itália tinha algumas observações de um asteroide, nomeado Ceres, o qual não foi visto depois de algumas semanas. O problema despertou o interesse de Gauss, pois a partir de um determinado número de observações queria descobrir a órbita do planeta. Assim, descobriria para onde os observadores deveriam apontar o telescópio. Para determinar a órbita de planetas, Gauss criou o método de mínimos quadrados que levou a resultados surpreendentes. O asteroide Ceres foi visualizado depois em posição muito próxima a indicada pelos cálculos realizados por Gauss (BOYER, 2010; EVES, 2004).

O objetivo do trabalho foi comparar os métodos iterativos de Newton e Gauss-Newton no ajuste do modelo Stanford e Smith (1972) aos dados de mineralização de carbono de palha de aveia no solo.
%Figuras (Figures)

%\begin{figure}[!htp]
%\begin{center}
%\includegraphics[scale=1]{xxx}
%\caption{xxxxx}\label{fig:1}
%\end{center}
%\end{figure}

%Tabelas: fonte Times New Roman tamanho 12. A identificação aparece na parte superior precedida da palavra designativa (tabela), seguida de seu número de ocorrência no texto (algarismos arábicos), e o respectivo título. A indicação da fonte aparece na parte inferior mesmo que seja de produção do próprio autor.

\section{Material e métodos}

Os dados utilizados para o ajuste do modelo foram extraídos de Giacomini et al. (2008). O experimento foi realizado em laboratório. Um argissolo vermelho distrófico arênico, da camada de 0-10cm de uma área que vinha sendo manejada em sistema plantio direto foi avaliado. A coleta da aveia foi feita no estádio de maturação fisiológica, submetida à secagem ao ar e armazenada em lugar seco até o momento da incubação. Neste trabalho foi avaliado um tratamento com palha de aveia incorporada no solo. A mineralização do C foi avaliada por meio da emissão de CO$_2$, durante a incubação, medindo-se a porcentagem de carbono mineralizado sempre nas mesmas unidades experimentais aos 3, 5, 9, 14, 20, 25, 30, 35, 45, 55, 65 e 80 dias do início da incubação.

Avaliou-se o modelo Stanford e Smith:
 $$C_i=C_0(1-exp(-kt_i))+\varepsilon_i$$
No modelo, $C_i$ corresponde a porcentagem do carbono adicionado mineralizado até o tempo $t_i$ (em dias), $i=1,2,...,n$; $C_0$ é o carbono potencialmente mineralizável; $k$ é a taxa de mineralização e $\varepsilon_i$ é o erro experimental (SILVA, et al., 2019a; SILVA, et al., 2019b).

O método de mínimos quadrados protosto por Gauss consiste em minimizar a soma de quadrados do erro. Para os modelos de regressão o erro é obtido por: $\mat{\varepsilon} = {\bf y}-f({\bf X},\mat{\beta})$ e a soma de quadrados do erro é dada por: 
\begin{align*}
\mat{\varepsilon}'\mat{\varepsilon} &= ({\bf y}-f({\bf X},\mat{\beta}))'({\bf y}-f({\bf X},\mat{\beta}))\\
&={\bf y}'{\bf y}-2f'({\bf X},\mat{\beta}){\bf y}+f'({\bf X},\mat{\beta})f({\bf X},\mat{\beta}).
\end{align*}

Observa-se que o ponto que minimiza a soma de quadrados do erro $\mat{\varepsilon}'\mat{\varepsilon}$ é obtido derivando-se a função em relação a $\mat{\beta}$ e igualando a zero, chega-se ao sistema de equações normais: $${\bf X}'f({\bf X},\mat{\beta})={\bf X}'{\bf y}$$
em que {\tiny ${\bf X}=\dfrac{\partial f({\bf X},\mat{\beta})}{\partial\mat{\beta}}$}. Logo, ${\bf X}$ é a matriz de derivadas parciais do modelo, e para modelo não linear a matriz ${\bf X}$ tem parâmetro pela definção de que um modelo é classificado como não linear se pelo menos uma derivada parcial da função depender de parâmetro. Portanto, $f({\bf X},\mat{\beta})$ e ${\bf X}$ dependem de $\mat{\beta}$. Deste modo, não é possível obter forma fechada para a solução $\hat{\mat{\beta}}$, sendo necessário o uso de método iterativo para aproximação da solução (DRAPER e SMITH, 2014).

Assim, precisa-se determinar ${\bf J}$ e ${\bf H}$ da função $\mat{\varepsilon}'\mat{\varepsilon}$ para usar o método de Newton para ter uma solução aproximada dos $\mat{\beta}$ que minimizam a soma de quadrados do erro:
$${\bf J}=\dfrac{\partial\mat{\varepsilon}'\mat{\varepsilon}}{\partial\mat{\beta}} = -2\dfrac{\partial f'({\bf X},\mat{\beta})}{\partial\mat{\beta}}{\bf y}+2\dfrac{\partial f'({\bf X},\mat{\beta})}{\partial\mat{\beta}}f({\bf X},\mat{\beta})=-2\dfrac{\partial f'({\bf X},\mat{\beta})}{\partial\mat{\beta}}({\bf y}-f({\bf X},\mat{\beta}))=-2{\bf X}'\mat{\varepsilon}$$
e
$${\bf H}=\dfrac{\partial^2\mat{\varepsilon}'\mat{\varepsilon}}{\partial\mat{\beta}\partial\mat{\beta}'} =\dfrac{\partial {\bf J}}{\partial\mat{\beta}}=-2\left[ {\bf X}'{\bf X} -\displaystyle\sum_{i=1}^{n}\dfrac{\partial^2f({\bf X}_i,\mat{\beta})}{\partial\mat{\beta}\partial\mat{\beta}'}\mat{\varepsilon_i} \right]$$

Pelo método de Newton:

$$\mat{\beta}^{n+1} = \mat{\beta^n} - ({\bf H}^n)^{-1}{\bf J}^n$$

Substituindo-se ${\bf J}$ e ${\bf H}$ na expressão acima, tem-se:

$$\mat{\beta}^{n+1} = \mat{\beta^n} -\left(-\dfrac{1}{2}\right)\left[ {\bf X}'{\bf X} -\displaystyle\sum_{i=1}^{n}\dfrac{\partial^2f({\bf X}_i,\mat{\beta})}{\partial\mat{\beta}\partial\mat{\beta}'}\mat{\varepsilon_i} \right]^{-1}(-2{\bf X}'\mat{\varepsilon})$$
ou
$$\mat{\beta}^{n+1} = \mat{\beta^n} - \left[ {\bf X}'{\bf X} -\displaystyle\sum_{i=1}^{n}\dfrac{\partial^2f({\bf X}_i,\mat{\beta})}{\partial\mat{\beta}\partial\mat{\beta}'}\mat{\varepsilon_i} \right]^{-1}{\bf X}'\mat{\varepsilon}$$

Para o seu problema de minimizar a soma de quadrados do erro $\mat{\varepsilon}'\mat{\varepsilon}$, Gauss desconsiderou o segundo termo da matriz Hessiana, pelo fato do termo {\tiny $\displaystyle\sum_{i=1}^{n}\dfrac{\partial^2f({\bf X}_i,\mat{\beta})}{\partial\mat{\beta}\partial\mat{\beta}'}\mat{\varepsilon_i}$} ser muito pequeno em relação a ${\bf X}'{\bf X}$ a medida que se aproxima da solução do sistema. Com essa otimização, criou o denominado método de Gauss-Newton que ficou da seguinte forma:

$$\mat{\beta}^{n+1} = \mat{\beta^n} - \left[ {\bf X}'{\bf X} \right]^{-1}{\bf X}'\mat{\varepsilon}$$
(BARD, 1974; JUDGE, et al., 1985; NOCEDAL e WRIGHT, 2006).

O método de Gauss-Newton também é denominado ``método da linearização", pois desconsiderar o termo do erro na matriz Hessiana equivale a considerar a equação $f({\bf X},\mat{\beta})$ a série de Taylor até o termo de primeira ordem (aproximação linear) e aplicar o método de mínimos quadrados (DRAPER e SMITH, 2014; MAZUCHELI e ACHCAR, 2002).

Considerou-se que o modelo ajustado convergiu se a diferença da soma de quadrados do erro (SQE) das iterações $i-1$ e $i$ foi menor que $10^{-4}$. Além disso, foi contado o número de iterações que o método precisou para convergir e calculou-se a SQE. Os valores iniciais para os parâmetros foram $\mat{\beta^0}' =(C_0^0 \quad k^0)'=(70 \quad 0,03)'$ As análise foram realizadas utilizando-se o software R (R Core Team, 2017). 

\section{Resultados e discussão}

Apresenta-se abaixo, com base no modelo Stanford e Smith, a matriz ${\bf X}$ usada na implementação dos algoritmos de Newton e Gauss-Newton e a matriz {\tiny $\dfrac{\partial^2f({\bf X}_i,\mat{\beta})}{\partial\mat{\beta}\partial\mat{\beta}'}$} usada na implementação do método iterativo de Newton. Uma das vantagens do método iterativo de Gauss-Newton é não precisar calcular a matriz {\tiny $\dfrac{\partial^2f({\bf X}_i,\mat{\beta})}{\partial\mat{\beta}\partial\mat{\beta}'}$} de segundas derivadas (NOCEDAL e WRIGHT, 2006).

\begin{center}
${\bf X}=  \begin{bmatrix} 
  1-e^{-kt_1} &  t_1C_0e^{-kt_1} \\ 
  \vdots  & \vdots \\ 
  1-e^{-kt_n} &  t_nC_0e^{-kt_n} \\
  \end{bmatrix}
$
e
$
\dfrac{\partial^2f({\bf X}_i,\mat{\beta})}{\partial\mat{\beta}\partial\mat{\beta}'} =  \begin{bmatrix} 
  0  & t_ie^{-kt_i} \\ 
  t_ie^{-kt_i} & -t_i^2C_0e^{-kt_i} \\
  \end{bmatrix}
$
\end{center}

Na Tabela 1 são apresentadas as estimativas dos parâmetros do modelo Stanford e Smith aos dados de mineralização do carbono de palha de aveia no solo considerando os métodos iterativos de Newton e Gauss-Newton. Observa-se que a diferença na estimativa do parâmetro $C_0$ foi na terceira casa decimal e do parâmetro $k$ na sexta casa decimal, o que na prática é desprezível, como ilustrado na Figura 1. Logo, ambos os métodos iterativos foram eficientes na estimatição dos parâmetros do modelo. Segundo Mazucheli e Achcar (2002) a escolha de bons valores iniciais são fundamentais para a convergência de qualquer método iterativo.\\

{\bf Tabela 1:} Estimativas dos parâmetros do modelo Stanford e Smith considerando os métodos iterativos (Fonte: próprio autor).
\begin{center}
\begin{tabular}{ccc}
\hline
Parâmetro & Método de Newton & Método Gauss-Newton \\
\hline
C$_0$ & 67,70202200 & 67,70131516 \\
k & 0,02414186 & 0,02414234 \\
\hline
\end{tabular}
\end{center}



\begin{center}
\includegraphics[scale=0.4]{figura}
\end{center}
{\bf Figura 1:} Ajuste do modelo Stanford e Smith considerando os métodos iterativos (Fonte: próprio autor).\\


Na Tabela 2 são apresentadas a análise da convergência considerando os métodos iterativos de Newton e Gauss-Newton, que precisou de apenas 6 e 4 iterações até atingir a convergência (SQE$_{i-1}$ $-$ SQE$_{i}$$<$$10^{-4}$), respectivamente, e a SQE apresentou diferença depois da quinta casa decimal (Tabela 2). O método Gauss-Newton é específico para estimar parâmetros de modelos de regressão não linear (MAZUCHELI e ACHCAR, 2002) e pela eficiência está implementado nos softwares até hoje (R Core Team, 2017).\\


{\bf Tabela 2:} Análise da convergência considerando os métodos iterativos (Fonte: próprio autor).
\begin{center}
\begin{tabular}{ccc}
\hline
%Parâmetro & Média {\it posteriori} & \multicolumn{2}{c}{HPD 95\%} \\
 & Método de Newton & Método Gauss-Newton \\
\hline
Número de iterações & 6 & 4\\
SQE & 27,35236 & 27,35236 \\
SQE$_{i-1}$-SQE$_{i}$ & 0,00006 & 0,00001 \\
\hline
\end{tabular}
\end{center}
SQE - soma de quadrados dos erros.\\



Para ilustrar os passos dos algoritmos, a Figura 2 foi feita com 6 iterações de cada método numérico. Percebe-se que a partir da quarta iteração o tamanho do passo foi muito pequeno.

\begin{center}
\includegraphics[scale=0.5]{traco}
\end{center}
{\bf Figura 2:} Passos dos métodos iterativos no ajuste do modelo Stanford e Smith aos dados de mineralização de palha de aveia (Fonte: próprio autor).\\

\newpage
Apresenta-se abaixo, com base no modelo Stanford e Smith, a matriz ${\bf X'X}$ e a matriz {\tiny $\displaystyle\sum_{i=1}^{n}\dfrac{\partial^2f({\bf X}_i,\mat{\beta})}{\partial\mat{\beta}\partial\mat{\beta}'}\mat{\varepsilon_i}$} na última iteração do método de Newton. Observa-se que os termos da matriz {\tiny $\displaystyle\sum_{i=1}^{n}\dfrac{\partial^2f({\bf X}_i,\mat{\beta})}{\partial\mat{\beta}\partial\mat{\beta}'}\mat{\varepsilon_i}$} são irrelevantes em relação aos termos da matrix ${\bf X'X}$ (NOCEDAL e WRIGHT, 2006).

\begin{center}
${\bf X'X}=  \begin{bmatrix} 
  3,4183 & 4967,5200  \\ 
  4967,5200 & 7869902,7300  \\
  \end{bmatrix}
$
e
$
\displaystyle\sum_{i=1}^{n}\dfrac{\partial^2f({\bf X}_i,\mat{\beta})}{\partial\mat{\beta}\partial\mat{\beta}'}\mat{\varepsilon_i} =  \begin{bmatrix} 
  0 & 0,0017 \\ 
 0,0017  & 47857,4700 \\
  \end{bmatrix}
$
\end{center}


\section{Conclusão}

O método dos mínimos quadrados proposto por Gauss tem grande aplicação na Estatística na estimação de parâmetros. O método Gauss-Newton, específico para estimar parâmetros de modelos não lineares, foi mais  eficiêntes que o método de Newton na estimação dos parâmetros do modelo Stanford e Smith ajustado a dados de mineralização de carbono. O método Gauss-Newton é implementado nos softwares para estimação de parâmetros de modelos de regressão não linear.

\section*{Agradecimentos}

Os autores agradecem a Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES) e ao Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq) pelo apoio nesta pesquisa.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                              %
%                  REFERENCES  (REFERÊNCIAS)                   %
%                                                              %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Referencias Bibliográficas}

\begin{flushleft}
\noindent BARD, Y. {\it Nonlinear parameter estimation}. 1a ed. Academic Press, 1974.\newline

\noindent BOYER, C. B. {\it História da Matemática,} 3 edição. São Paulo: Bl\"{u}cher, 2010.
\newline

\noindent DRAPER, N. R.; SMITH, H. {\it Aplied regression analisys}. 3a ed, reprint. John Wiley, 2014.\newline

\noindent EVES, H. W. {\it Introdução a história da matemática.} São Paulo: UNICAMP, 2004.
\newline

\noindent GIACOMINI, S. J.; AITA, C.; MIOLA, E. C. C.; RECOUS, S. Mineralização do carbono da palha de aveia e dejetos de suínos aplicados na superfície ou incorporados ao solo. {\it Revista Brasileira de Ciência do Solo,} v. 32, p. 2661-2668, 2008.\newline

\noindent JUDGE, G. G.; GRIFFITHS, W. E.; HILL, R. C.; L\"{U}TKEPOHL, H.; LEE, T. {\it The theory and practice of econometrics}. 2a ed. John Wiley, 1985.\newline

\noindent MAZUCHELI, J.; ACHCAR, J. A. Algumas considerações em regressão não-linear. {\it Acta Scientiarum,} v. 24, p. 1761-1770, 2002.\newline

%\noindent MACHADO, E. J. et al. Estimação de um modelo de espécies de macroinvertebrados bentônicos via análise bayesiana do modelo de Michaelis-Menten. {\it Revista Brasileira de Biometria,} v. 30, p. 106-123, 2012.\newline

%\noindent PEREIRA, J. M. et al. Comparação entre modelos para predição do nitrogênio mineralizado: uma abordagem bayesiana, {\it Ciência e Agrotecnologia,} v. 33, p. 1792-1797, 2009.\newline

\noindent NOCEDAL, J. WRIGHT, S. {\it Numerical Optimization}. 2a ed. Springer, 2006.\newline


\noindent R CORE TEAM. {\it R: A language and environment for statistical computing.} R Foundation for Statistical Computing, Vienna, Austria. 2017. URL \url{http://www.R-project.org/}.\newline

\noindent SILVA, E. M.; RIBEIRO, T. D.; FERNANDES, J. G.; MUNIZ, J. A. Descrição da mineralização do carbono de dejetos de suínos e palha de aveia no solo por modelos não lineares. {\it Revista Agrogeoambiental,} v. 11, p. (no prelo), 2019a.\newline

\noindent SILVA, E. M.; SILVEIRA, S. C.; RIBEIRO, T. D.; MUNIZ, J. A. Ajuste da decomposição do lodo de esgoto e palha de aveia por modelos não lineares. {\it Revista Agrogeoambiental,} v. 11, p. (no prelo), 2019b.\newline

%\noindent SOUZA, E. M. Non-linear modeling of zinc extracted from a sewage sludge-treated. {\it Acta Scientiarum,} v. 32, p. 193-199, 2010.\newline

\noindent STANFORD, G.; SMITH, S, J. Nitrogen mineralization potentials of soil. {\it Soil Science Society of America Jornal,} v. 36, p. 465-471, 1972.\newline

%\noindent ZEVIANI, et al. Modelos não lineares para a liberação de potássio de estercos animais em latossolos. {\it Ciência Rural,} v. 42, p. 1789-1796, 2012.\newline

\end{flushleft}

\end{document}

