Estimativa da relação hipsométrica em clones de Eucalyptus sp. com o modelo de curtis ajustado por métodos bayesianos empíricos

Estimate of the hypsometric relationship for Eucalyptus sp. clones with the curtis' model fitted by empirical bayesian methods

Estimativa da relação hipsométrica em clones de Eucalyptus sp. com o modelo de curtis ajustado por métodos bayesianos empíricos

Estimate of the hypsometric relationship for Eucalyptus sp. clones with the curtis' model fitted by empirical bayesian methods

Cláudio Roberto Thiersch^I; Marinho Gomes de Andrade^II; Monica Fabiana Bento Moreira^III; Selene Loibel^IV

^IUniversidade Federal de São Carlos , DCA. E-mail: <[e-mail]>

^IIUniversidade de São Paulo, ICMC-SME. E-mail: <[e-mail]>

^IIIUniversidade Federal de São Carlos , Cac-Sor. E-mail: <[e-mail]>

^IVUniversidade Estadual de São Paulo, DEMAC. E-mail: <[e-mail]>

RESUMO

Neste trabalho foi considerado o modelo de Curtis para a relação hipsométrica em clones de Eucalyptus sp. com os parâmetros sujeitos a restrições. Para fazer a inferência dos parâmetros do modelo com restrições, utilizou-se uma abordagem bayesiana com densidade a priori construída empiricamente. As estimativas bayesianas são calculadas com a técnica de simulação de Monte Carlo em Cadeia de Markov (MCMC). O método proposto foi aplicado a diferentes conjuntos de dados reais, dos quais foram selecionados cinco para exemplificar os resultados. Estes foram comparados com os resultados obtidos pelo método de mínimos quadrados, destacando-se a superioridade da abordagem bayesiana proposta.

Palavras-chave: Modelo de regressão não linear, Inferência bayesiana, Restrições nos parâmetros, MCMC.

ABSTRACT

The model of Curtis was considered in this work for the hypsometric relationship in Eucalyptussp. clones with parameters submitted to restrictions. A Bayesian a priori density was empirically constructed to infer the parameters of the models with restrictions. The Bayesian estimates are calculated using Monte Carlo Markov Chain (MCMC) Method. The proposed method was applied to different groups of real data from which five were selected to show the results. Those were compared to the achieved results by the minimum square method and the superiority of the Bayesian approach is highlighted.

Keywords: Curtis' Model, Bayesian inference, Height-diameter relationship, Restriction parameters, MCMC.

1. INTRODUÇÃO

No inventário florestal, a quantificação das alturas das árvores é fundamental, principalmente, para a estimativa do volume de madeira dos povoamentos e para fazer classificações destes quanto às suas produtividades. Por sua vez, a altura é uma medida indireta sujeita a erros e, normalmente, implica maior tempo para sua quantificação quando comparada com o diâmetro, fato que torna a operação onerosa.

Para contornar esses problemas, tem-se como prática medir na parcela a altura de parte das árvores e os diâmetros de todas estas, bem como através dos pares altura-diâmetro mensurados, estabelecer relação matemática que possibilite estimar as alturas das demais árvores contidas nas parcelas. Essas relações matemáticas são conhecidas como equações hipsométricas, as quais buscam estimar as alturas das árvores através da relação Dap (diâmetro do tronco da árvore a 1,3 m do solo) e H (altura total da árvore).

A grande dificuldade da escolha do melhor modelo para representar essas relações hipsométricas se deve à não linearidade da relação entre as variáveis envolvidas e as restrições impostas aos parâmetros dos modelos, por razões biológicas (PLÁCIDO et al., 2004; SOARES et al., 2004). Esses dois fatores aliados tornam o ajuste de modelos estatísticos para representar as relações hipsométricas um problema de regressão não linear, com restrições nos parâmetros.

Para que as estimativas obtidas com essas equações sejam precisas, os métodos e modelos adotados tornam-se cada vez mais sofisticados (MACHADO et al., 1993; SOARES et al., 2004). Por sua vez, a abordagem do problema com os métodos convencionais de inferência estatística, como método de máxima verossimilhança e método de mínimos quadrados, torna-se ineficiente para considerar as restrições nos parâmetros de vários modelos propostos.

Neste trabalho foi considerado um desses modelos propostos na literatura, conhecido como o modelo de Curtis (CURTIS, 1967). O modelo de Curtis é relativamente simples devido à possibilidade de linearização do modelo. Mas a presença de restrições nos parâmetros torna o problema de inferência complexo.

Para resolver o problema de inferência com restrições nos parâmetros para o modelo de Curtis, é proposta neste trabalho uma abordagem bayesiana. Os métodos bayesianos consistem basicamente em assumir que os parâmetros de um modelo são variáveis aleatórias e as informações a priori sobre esses parâmetros podem ser modeladas usando distribuições de probabilidade. As distribuições a priori são combinadas com as informações provenientes de uma amostra, dada pela função de verossimilhança, resultando em uma distribuição a posteriori para os parâmetros do modelo. A grande vantagem dessa abordagem consiste na facilidade com que certos problemas de inferência podem ser tratados, como é o caso das relações hipsométricas com restrições nos parâmetros.

2. MATERIAL E MÉTODOS

2.1. Localização e caracterização da área

Os dados utilizados para testar a metodologia proposta foram coletados em florestas plantadas de clones de Eucalyptus sp., da Votorantim Celulose e Papel (VCP), na região de Itapetininga, no Sul do Estado de São Paulo, com coordenadas geográficas 23º562 S de latitude e 42º212 W de longitude, a uma altitude média de 655 m.

De acordo com Golfari et al. (1978), nessa região a precipitação média anual é de 1.600 mm e a temperatura média anual, de 20 ºC. Os solos são predominantemente Latossolo Vermelho distrófico.

2.2. Amostragem

Foram utilizadas cerca de 1.000 parcelas locadas em povoamentos de Eucalyptus sp. com idades variando de 2 a 8 anos, das quais foram selecionadas cinco parcelas para exemplificar os resultados. A parcela A foi mensurada em um povoamento de 2,33 anos de idade e espaçamento de 3,0×3,0 m; a parcela B em povoamento de 3,46 anos e espaçamento de 4,0×2,2 m; a parcela C em povoamento de 4,94 anos e espaçamento de 5,0×1,5 m; a parcela D em povoamento de 6,12 anos e espaçamento de 3,0×3,0 m; e a parcela E em povoamento de 7,64 anos, com espaçamento de 5,0×1,8 m.

Nessas unidades amostrais, cujo formato era circular e com área de 500 m², foram medidos os diâmetros a 1,3 m do solo de todas as árvores e a altura total de 10 árvores centrais e de cinco dominantes.

2.3. Modelo de Curtis

O modelo de Curtis para representar a relação hipsométrica é dado por:

em que H_i representa a altura total da árvore i; Dap é o diâmetro do tronco da árvore i a uma altura de 1,30 m do solo; e b₀ e b₁ são os parâmetros do modelo que, por razões biológicas, estão sujeitos às restrições do tipo β₀ > 0 e β₁ > 0, sendo ε_i o erro associado do modelo.

2.4. Estimador de mínimos quadrados

Considerando o conjunto de observações D={(H_i,Dap_i), i=1, ...,n}, sob a hipótese de que a relação entre altura H e diâmetro Dap para essa amostra é dada pelo modelo de Curtis, pode-se escrever para cada par (H_i,Dap_i) a relação:

em que y_i é o ln, x_i é 1/Dap_i, e_i é ln(ε_i) e os parâmetros e já foram descritos anteriormente. Considerando as n observações, pode-se escrever a relação (2) na forma matricial como:

com Y = (y_i , ..., y_n)' e X = (1' ,-x'), sendo 1=(1,...,1)' um vetor (n×1) somente de uns e x = (x',...,x_n)'. Além disso, tem-se e = (e_i, ...,e_n)' e β = (β₀,β₁ )'. Portanto, o estimador de mínimos quadrados para o vetor de parâmetros â irrestrito pode ser calculado de forma direta por:

A matriz variância do erro e pode ser expressa como σ ²= E(e' e)/n. Portanto, um estimador não viesado para a variância σ ² pode ser calculado como:

Devido à ocorrência natural em determinados povoamentos de algumas árvores em que não há proporcionalidade direta entre os diâmetros a 1,3 m do solo e as alturas totais e, também, devido ao pequeno tamanho (n) das amostras que geralmente são utilizadas em situações práticas, as estimativas de mínimos quadrados irrestritas para os parâmetros dadas em (4) resultam, em muitos casos, em valores inconsistentes biologicamente, ou seja, estimativas negativas que não permitem predizer corretamente a relação entre altura e diâmetro das árvores (THIERSCH, 2007). Para resolver esse problema, a solução é a inclusão das restrições no problema de estimação dos parâmetros, e isso leva a considerar, juntamente com a equação (3), as restrições do tipo β ∈ R, em que R={a_i<β_i<b_i,i=0,...,p},p ∈ N. Na definição dessas restrições são considerados, inclusive, os casos em que alguns ou todos os a_i podem ser -∞ e os b_i, + ∞. A solução numérica desse problema pode ser encontrada, mas exige condições iniciais adequadas; caso contrário, implica muitas dificuldades de convergência.

Para contornar tal dificuldade, está sendo proposta uma abordagem bayesiana para o cálculo dessas estimativas. Na abordagem bayesiana, as informações a priori sobre os parâmetros podem compensar, parcialmente, a pouca informação extraída somente dos dados quando estes constituem uma pequena amostra.

2.5. Abordagem bayesiana do modelo de Curtis

A abordagem bayesiana do problema de inferência dos parâmetros de um modelo assume que esses parâmetros são variáveis aleatórias, e qualquer informação inicial sobre eles pode ser modelada por uma função densidade de probabilidade a priori. Combinando essas densidades a priori com a função de verossimilhança dos dados, por meio do Teorema de Bayes, chega-se à distribuição a posteriori conjunta. Denotando o vetor de parâmetros de um modelo por θ, a densidade a priori por π₀(θ) e a função de verossimilhança associada a um conjunto de observações D por L(D|θ), então a distribuição a posteriori é dada por:

A integral em (6) é uma integral múltipla sobre o domínio de definição dos parâmetros θ ∈ , que representa a constante normalizadora da distribuição a posteriori, portanto é uma função somente dos dados D observados. É comum adotar a notação de proporcionalidade para representar a densidade a posteriori, dada por:

Para maiores detalhes sobre inferência bayesiana, ver Box e Tiao (1973).

a) Densidade de probabilidade a priori empírica

Para considerar as restrições do tipo β ∈ R, em que R={a_i<β_i<b_i, i=1, ..., n}, serão levadas em conta as densidades a priori conjugadas Normal-Gama Truncada para os parâmetros assim β|τ∼|NT(β₀,(τP)^-1), em que τ = 1/σ²>τ∼G(ν₀ ,λ 0), ou seja, tem-se:

Em (8), µ = (µ₀, µ₁)é o vetor de hiperparâmetros de locação e P é a matriz de precisão, ambos conhecidos. I_R(β) é uma função indicadora tal que I_R(β)=1 se β ∈ R e I_R(β)=0, caso contrário. A densidade a priori para τ é dada por:

Em (9) ν,₀ e λ₀ são hiperparâmetros também conhecidos. A priori conjunta é dada por:

Neste trabalho foi adotada uma abordagem empírica para determinação dos hiperparâmetros µ e P. Para isso, consideraram-se os percentis (H_p,Dap_p) e (H_q,Dap_q), com p<q encontrados em dados históricos coletados anteriormente às observações da amostra D, os quais não ocorreram obrigatoriamente de forma pareada (neste trabalho foi considerado nas aplicações p=25% e q=75%, mas outros percentis podem ser considerados sem perda de generalidade). Assumindo que o modelo de Curtis (1) é válido para esses percentis, pode-se escrever:

em que estamos considerando y_k=ln(H_k),x_k=ln(Dap_k) para k=p,q e m = (µ₀^*, µ₁^*). A matriz de precisão P também pode ser encontrada de forma empírica, considerando-se um coeficiente de variação constante CV_j = γ_j|µ_j, j=0,1 (BROEMELING, 1985), que expresse a confiança nos hiperparâmetros µ. Dessa forma, tem-se:

Assume-se também que β ₀ e β ₁ são independentes a priori, de forma que P_i,j=0 para i ≠ j, (i,j=0,1).

b) A densidade de probabilidade a posteriori

Utilizando a metodologia bayesiana (ver detalhes em ANDRADE et al., 2007), encontra-se a densidade de probabilidade a posteriori para β e τ, que é dada por:

em que W=Y' Y+V, V=(X'X+P) e é o estimador bayesiano irrestrito, dado por:

sendo =(X' X)^-1 X' Y o estimador de mínimos quadrados.

Para encontrar os estimadores bayesianos a partir da densidade a posteriori conjunta (14), adotou-se um algoritmo de simulação de Monte Carlo em Cadeia de Markov (GILKS, 1998). Para isso, geraram-se amostras das densidades condicionais a posteriori, π(β|τ,D), que é Normal-Truncada e π(τ|β,D), que é uma Gama, ambas dadas por:

em que Λ(β)=(β -)' V(β -)+W. O algoritmo amostrador de Gibbs Sampler utilizado para gerar amostras da densidade a posteriori pode ser visto detalhadamente em Gilks (1998), aqui ele é apresentado de forma simplificada por:

• Algoritmo Gibbs Sampler:

1. Fazer j=0 e considerar a condição inicial, τ⁽⁰⁾.

2. Gerar β^(j+1) de β^j+1)|τ^(j),D∼ NT(^(j+1)V)^-1).

3. Gerar τ ^(j+1) de τ ^(j+1)|β^(j+1),D ∼ G(n/2+ ν₀ +1,Λ(β^(j+1))+2λ₀).

4. Fazer j←j+1 e repetir os passos (2), (3) e (4).

Para encontrar os estimadores bayesianos obtidos via simulação de Monte Carlo, adotou-se o procedimento sugerido em Gilks (1998). Considerou-se um período de aquecimento (burn in) de 50% das iterações iniciais. Os passos (2), (3) e (4) do algoritmo Gibbs Sampler são repetidos até que a convergência monitorada graficamente e com o critério de Geweke (1992) seja verificada. Então, uma amostra é selecionada a cada 10 iterações. Denotando por θ ^(j)={β₀^(j), β₁^(j), β_0,R^(j), β_1,R^(j), τ^(j)}, j=1,...,M a amostra resultante do algoritmo Gibbs Sampler, os estimadores bayesianos obtidos via simulação de Monte Carlo são dados por:

E os intervalos de credibilidade [_k^MC(α), _k^MC(1-(α)] são estimados com base nos percentis da amostra, dados por:

P[_k^MC(α)≤_k^MC(1-(α)]=1-(α)

3. RESULTADOS E DISCUSSÕES

3.1. Densidades a priori empíricas

Os parâmetros das densidades a priori empíricas consideradas para cada parcela são dados na Tabela 1.

Esses valores foram calculados usando-se as equações (11 a 13) com dados de parcelas já analisados.

Consideraram-se nas equações (11) e (12) os percentis de p=25% e q=75% e na equação (13), fizeram-se os coeficientes de variação CV₀=1 e CV₁=10. Esses valores são arbitrários, porém a escolha de coeficientes de variações grandes é sempre recomendável para garantir que as densidades a priori sejam curvas aplainadas (flat curves).

3.2. Análise a posteriori

As estimativas bayesianas foram obtidas via simulação de Monte Carlo em Cadeia de Markov (MCMC). O algoritmo Gibbs Sampling gerou uma amostra de tamanho 20.000, da qual foram descartadas as primeiras 10.000 (burn in). Selecionou-se da segunda parte um valor a cada 10 valores gerados (ou seja, 10.001, 10.010, 10020,..., 20.000), perfazendo, assim, uma amostra de tamanho 1.000, que foi utilizada para o cálculo das estimativas bayesianas via MCMC.

Um sumário das estimativas bayesianas e de mínimos quadrados encontra-se na Tabela 2. Comparando as estimativas bayesianas com as de mínimos quadrados sem restrições apresentadas na Tabela 2, notou-se que várias estimativas de mínimos quadrados apresentam valores negativos, violando, portanto, as restrições biológicas impostas aos parâmetros. No entanto, as estimativas bayesianas, dadas pela média das amostras geradas da densidade a posteriori via MCMC, apresentam sempre valores positivos. Os intervalos de credibilidade bayesiana IC(95%) apresentados na Tabela 2 revelam a precisão dessas estimativas. Nota-se na que, quando os estimadores de Mínimos Quadrados (MQ) resultam em valores positivos, essas estimativas coincidem com as estimativas bayesianas, como mostrado no resultado da parcela B.

As densidades marginais a posteriori obtidas via MCMC, juntamente com as densidades a priori empíricas, são apresentadas na Figura 1, para cada um dos parâmetroTabela 2s do modelo de Curtis ajustado aos dados da amostra da parcela E. Os resultados das outras parcelas são semelhantes.

Analisando as densidades apresentadas na Figura 1, observou-se que a combinação das informações fornecidas a priori com as informações extraídas dos dados através da função de verossimilhança localiza o domínio de definição das densidades a posteriori dentro do domínio das densidades a priori empíricas, definido pelas restrições impostas de forma biologicamente aceitável. Além disso, a análise a posteriori aumentou a precisão das estimativas dos parâmetros, e isso pode ser constatado nos gráficos da Figura 1, observando-se densidades a posteriori mais fechadas em relação às densidades a priori. Esses dois aspectos (localização e forma das densidades a posteriori) destacaram a validade da proposta empírica de construção das densidades a priori e também a vantagem da abordagem bayesiana para tratar o problema de inferência dos parâmetros e das relações hipsométricas com restrições nos parâmetros.

As curvas estimadas usando os estimadores bayesianos e de mínimos quadrados para cada uma das amostras dadas na Tabela 2 são apresentadas na Figura 2. A análise dessas curvas revela que, em conformidade com as estimativas apresentadas na Tabela 2, as estimativas de mínimos quadrados (MQ) resultam em uma relação biologicamente inconsistente para representar a relação hipsométrica nas amostras das parcelas A, C, D e E. Na amostra da parcela B, ambas as abordagens coincidem, mostrando que há consistência da abordagem bayesiana proposta com a de mínimos quadrados, quando esta última é factível.

4. CONCLUSÕES

A abordagem bayesiana empírica proposta em conjunto com a técnica de solução via MCMC mostrouse vantajosa na consideração das restrições dos parâmetros do modelo. As estimativas obtidas pelo método de mínimos quadrados (MQ) resultam em uma relação biologicamente inconsistente para representar a relação hipsométrica em quatro dos cinco conjuntos de dados analisados, enquanto as estimativas bayesianas respeitaram os limites impostos aos parâmetros para todos os conjuntos de dados analisados.

A modelagem bayesiana proposta é, portanto, recomendável como técnica alternativa precisa e eficiente para a estimação dos parâmetros das relações hipsométricas. No entanto, devido ao esforço computacional dos métodos MCMC, uma avaliação da abordagem proposta neste artigo para tratar uma grande massa de dados, como geralmente ocorre na prática, é uma questão a ser considerada cuidadosamente.

5. REFERÊNCIAS

Recebido em 04.09.2007 aceito para publicação em 19.12.2012.

• ANDRADE, M. G. et al. Uma abordagem bayesiana empírica para o modelo hipsométrico de curtis com restrições nos parâmetros: uma aplicação para clones de Eucalyptus sp. São Carlos: USP, 2007. p.15. (Relatório Técnico ICMC).
• BOX, G. E.; TIAO, G. C. Bayesian inference in statistical analysis New York: Addison-Wesley, 1973.
• BROEMELING, L. D. Bayesian analysis of linear models New York: Marcel Dekker, 1985. 545p.
• CURTIS, R. O. Height diameter and height diameter age equations for second growth Douglas-fir. Forest Science, v.13, n.4, p.365-375, 1967.
• GEWEKE, J. Evaluating the accuracy of samplingbased approaches to calculating posterior moments. In: BERNARDO, J. O. et al. (Ed.) Bayesian statistics 4 Oxford: Clarendon Press, 1992.
• GILKS, W. R.; RICHARDSON, S.; SPIEGELHALTER, D. J. Markov Chain Monte Carlo in practice New York: Chapman & Hall, 1998. 486p.
• GOLFARI, L.; CASER, R. L.; MOURAA, V. P. G. Zoneamento ecológico esquemático para reflorestamento no Brasil: 2Ş aproximação. Belo Horizonte: Centro de Pesquisas Florestais da Região do Cerrado, 1978. 66p.
• MACHADO, S. A.; BASSO, S. F.; BEVILACQUA JR, V. G. Teste de modelos matemáticos para o ajuste da relação hipsométrica em diferentes sítios e idades para plantações de Pinus elliottii no Estado do Paraná. In: CONGRESSO FLORESTAL PANAMERICANO, 1.; CONGRESSO FLORESTAL BRASILEIRO, 7.,Curitiba, 1993. Anais..., São Paulo: Sociedade Brasileira de Silvicultura, 1993. v.2. p.553-556.
• PLÁCIDO, A. C. et al. Dinâmica da relação hipsométrica em função da idade, do sítio e da densidade inicial de povoamentos de Bracatinga da Região Metropolitana de Curitiba, PR. Revista Árvore, v.28, n.4, p.517-533, 2004.
• SOARES, T. S. et al. Uso de diferentes alternativas para viabilizar a relação hipsométrica no povoamento florestal. Revista Árvore, v.28, n.6, p.845-854, 2004.
• THIERSCH, C. R. Uma nova proposta metodológica para cubagem rigorosa e ajuste de modelos hipsométricos 2007. 75f. Tese (Doutorado em Engenharia Florestal) - Universidade Federal de Lavras, Lavras, 2007.

Datas de Publicação

• Publicação nesta coleção
  28 Mar 2013
• Data do Fascículo
  Fev 2013

Histórico

• Recebido
  04 Set 2007
• Aceito
  19 Dez 2012