%\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}