( PDF ) Comparação entre modelos para predição do nitrogênio mineralizado: uma abordagem bayesiana

Download PDF

ads:

COMPARAÇÃO ENTRE MODELOS PARA

PREDIÇÃO DO NITROGÊNIO

MINERALIZADO: UMA ABORDAGEM

BAYESIANA

JANSER MOURA PEREIRA

2006

ads:

Livros Grátis

http://www.livrosgratis.com.br

Milhares de livros grátis para download.

JANSER MOURA PEREIRA

COMPARAÇÃO ENTRE MODELOS PARA PREDIÇÃO DO

NITROGÊNIO MINERALIZADO: UMA ABORDAGEM BAYESIANA

Tese apresentada à Universidade Federal de Lavras

como parte das exigências do Programa de Pós-

Graduação em Agronomia, área de concentração em

Estatística e Experimentação Agropecuária, para a

obtenção do título de “Doutor”.

Orientador

Prof. Dr. Joel Augusto Muniz

Co-orientadora:

Profa Dra Thelma Sáfadi

LAVRAS

MINAS GERAIS – BRASIL

2006

ads:

Ficha Catalográfica Preparada pela Divisão de Processos Técnicos da

Biblioteca Central da UFLA

Pereira, Janser Moura

Comparação entre modelos para predição do nitrogênio

mineralizado: uma abordagem bayesiana / Janser Moura Pereira.

–

Lavras: UFLA, 2006.

79 p. : il.

Orientador: Joel Augusto Muniz.

Tese (Doutorado) – UFLA.

Bibliografia.

1. Inferência bayesiana. 2. Modelos não lineares. 3. Amostrador

de Gibbs. 4. Metropolis Hastings. 5. Fator de Bayes. 6. Critério de

Informação Bayesiano.

I. Universidade Federal de Lavras. II.

Título.

CDD-519.542

JANSER MOURA PEREIRA

COMPARAÇÃO ENTRE MODELOS PARA PREDIÇÃO DO

NITROGÊNIO MINERALIZADO: UMA ABORDAGEM BAYESIANA

Tese apresentada à Universidade Federal de Lavras

como parte das exigências do Programa de Pós-

Graduação em Agronomia, área de concentração em

Estatística e Experimentação Agropecuária, para a

obtenção do título de “Doutor”.

APROVADA em 15 de dezembro de 2006.

Prof. Dr. Marcelo Tavares UFU

Profa Dra Thelma Sáfadi UFLA

Prof. Dr. Carlos Alberto Silva UFLA

Prof. Dr. Renato Ribeiro de Lima UFLA

Prof. Dr. Joel Augusto Muniz

UFLA

(Orientador)

LAVRAS

MINAS GERAIS – BRASIL

Aos meus pais, Alaor Alves de Moura e Cleuza Helena Pereira, pelo amor e

incentivo.

À minha esposa, Andréa Clélia, pelo seu carinho, amizade e

amor em todos os momentos.

Aos meus irmãos, Jander Antônio e Tatiana Anselma,

pela amizade e carinho.

DEDICO

AGRADECIMENTOS

A DEUS, por tudo.

À Universidade Federal de Lavras (UFLA), ao Programa de Pós-

Graduação em Estatística e Experimentação Agropecuária e ao Departamento de

Ciências Exatas (DEX), pela oportunidade e confiança no meu trabalho.

Ao professor Joel Augusto Muniz e à professora Thelma Sáfadi, pela

orientação, atenção e amizade.

Aos professores membros da banca examinadora, pelas sugestões e

colaboração.

Aos professores do Departamento de Ciências Exatas da UFLA, pela

amizade e pelos ensinamentos.

Aos funcionários do Departamento de Ciências Exatas da UFLA, pela

assistência e eficiência com que nos atenderam.

Aos colegas de curso pós-graduação, pelo incentivo, companheirismo e

fraternidade, em especial aos amigos Denismar e José Waldemar.

À Coordenação de Aperfeiçoamento de Pessoal de Nível Superior

(Capes), pela concessão da bolsa de estudos.

A todos que, direta ou indiretamente, contribuíram para a realização

deste trabalho.

SUMÁRIO

Página

LISTA DE TABELAS............................................................................ i

LISTA DE FIGURAS............................................................................. ii

RESUMO................................................................................................ iii

ABSTRACT............................................................................................ iv

1 INTRODUÇÃO.................................................................................. 1

2 REFERENCIAL TEÓRICO.............................................................. 3

2.1 Comportamento do fenômeno: mineralização do nitrogênio........... 3

2.2 Modelos não lineares........................................................................ 4

2.2.1 Modelos não lineares adotados...................................................... 6

2.3 Inferência bayesiana.......................................................................... 8

2.3.1 Obtenção das distribuições posterioris condicionais completas 11

2.3.2 Priori não informativa.................................................................... 11

2.3.3 Intervalo de credibilidade............................................................... 12

2.3.3.1 Intervalo de máxima densidade a posteriori ou intervalo

HPD.........................................................................................

2.3.4 Erro de Monte Carlo...................................................................... 13

3 MATERIAL E MÉTODOS................................................................. 15

3.1 Material............................................................................................. 15

3.2 Métodos............................................................................................. 15

3.2.1 Métodos de reamostragem............................................................. 16

3.2.1.1 Método da rejeição...................................................................... 16

3.2.1.2 Método de reamostragem ponderada.......................................... 18

3.2.2 Método de Monte Carlo via cadeias de Markov (MCMC)............ 19

3.2.2.1 Amostrador de Gibbs.................................................................. 19

3.2.2.2 Algoritmo de Metropolis-Hastings............................................. 20

3.3 Diagnósticos de convergência........................................................... 21

3.3.1 Critério de Raftery & Lewis........................................................... 22

3.3.2 Critério de Geweke........................................................................ 23

3.4 Comparação dos modelos não lineares............................................. 24

3.4.1 Fator de Bayes (FB)....................................................................... 24

3.4.2 Critério de Informação Bayesiano (BIC)....................................... 25

3.5 Posteriori conjunta............................................................................ 25

3.5.1 Posteriori conjunta para o modelo não linear de Stanford &

Smith...........................................................................................

3.5.2 Posteriori conjunta para o modelo não linear de Marion............... 26

3.5.3 Posteriori conjunta para o modelo não linear de Cabrera.............. 28

4 RESULTADOS E DISCUSSÃO......................................................... 30

4.1 Obtenção das posterioris condicionais completas para o modelo de

Stanford & Smith...........................................................................

4.2 Obtenção das posterioris condicionais completas para o modelo de

Marion............................................................................................

4.3 Obtenção das posterioris condicionais completas para o modelo de

Cabrera...........................................................................................

4.4 Análise dos dados.............................................................................. 39

5 CONCLUSÕES................................................................................... 49

6 REFERÊNCIAS BIBLIOGRÁFICAS................................................. 51

7 ANEXO................................................................................................ 56

LISTA DE TABELAS

TABELA

Página

TABELA 1

Fator de dependência do critério de Raftery &

Lewis (FD) e p-valor do critério de Geweke dos

modelos (2.1), (2.2) e (2.3).......................................

TABELA 2

Média e desvio padrão a posteriori dos parâmetros

dos modelos (2.1), (2.2) e (2.3), com seus

respectivos intervalos HPD 95%, erro de Monte

Carlo (Erro MC).......................................................

TABELA 3

Fator de Bayes (FB) dos modelos (2.1), (2.2) e

(2.3)........................................................................... 44

TABELA 4

Critério de Informação Bayesiano (BIC) dos

modelos (2.1), (2.2) e (2.3)....................................... 44

LISTA DE FIGURAS

FIGURA

Página

FIGURA 1

N potencialmente mineralizável (

), com base

nas quantidades acumuladas de N mineralizado

durante a incubação sucessiva..................................

FIGURA 2

Trajetórias e traços das cadeias dos parâmetros

do modelo de Stanford & Smith.............

FIGURA 3

Trajetórias e traços das cadeias dos parâmetros

Nkb

do modelo de Marion..........................

FIGURA 4

Trajetórias e traços das cadeias dos parâmetros

Nk k

11 0

do modelo de Cabrera.......................

FIGURA 5

Densidades e histogramas dos parâmetros

do modelo de Stanford & Smith.............

FIGURA 6

Densidades e histogramas dos parâmetros

Nkb

do modelo de Marion..........................

FIGURA 7

Densidades e histogramas dos parâmetros

Nk k

11 0

do modelo de Cabrera.......................

RESUMO

PEREIRA, Janser Moura. Comparação entre modelos para predição do

nitrogênio mineralizado: uma abordagem bayesiana. 2006. 79 p. Tese

(Doutorado em Estatística e experimentação Agropecuária) – Universidade

Federal de Lavras, Lavras, MG.

Estudos recentes utilizam a inferência bayesiana nas mais diversas áreas.

Portanto, propõe-se neste trabalho desenvolver uma abordagem bayesiana para

predizer as quantidades de nitrogênio mineralizados por intermédio de modelos

não lineares, ou seja, ajustar um modelo de probabilidade para um grupo de

dados e sintetizar o resultado por meio de uma distribuição de probabilidade

para os parâmetros dos modelos. Os modelos não lineares considerados para

avaliar a dinâmica da mineralização do nitrogênio e para ilustrar o procedimento

bayesiano foram: modelo de Stanford & Smith, modelo de Marion e modelo de

Cabrera. A comparação dos modelos ocorreu por meio do Fator de Bayes (FB) e

do Critério de Informação Bayesiano (BIC). Neste trabalho, foram utilizados o

amostrador de Gibbs e o Metropolis Hastings, para inferência dos parâmetros.

Implementou-se um algoritmo no software R para realizar a amostragem de

Gibbs das distribuições posterioris dos parâmetros dos modelos. A convergência

das cadeias foram monitoradas por meio de análise gráfica e pelos critérios de

Geweke e Raftery & Lewis que estão implementados no pacote BOA,

executável no software R. O modelo de Cabrera foi o que proporcionou melhor

qualidade de ajuste ao conjunto de dados de mineralização de nitrogênio,

seguido pelo modelo de Stanford & Smith e, por último, o de Marion. Em

virtude dos resultados apresentados, pode-se atestar que a metodologia

bayesiana apresentou bons resultados na estimação dos parâmetros dos modelos,

ou seja, o ajuste de modelos por meio de distribuições posteriori condicionais

completas constitui uma metodologia confiável e bastante precisa.

Palavras Chaves: inferência bayesiana, modelos não lineares, amostrador de

Gibbs, Metropolis Hastings, fator de Bayes, critério de informação bayesiano.

Comitê Orientador: Prof. Dr. Joel Augusto Muniz – UFLA (Orientador) e

Profa Dra Thelma Sáfadi – UFLA (Co-orientadora)

iii

ABSTRACT

PEREIRA, Janser Moura. Comparison of the models for prediction of the

mineralized nitrogen: a bayesian boarding. 2006. 79 p. Thesis (Doctorate in

Statistics and Agricultural Experimentation) – Federal University of Lavras,

Lavras, Minas Gerais, Brazil.

Recent studies use the Bayesian Inference in the most several areas.

Therefore, intends in this work to develop a Bayesian boarding to predict of the

nitrogen mineralized through nonlinear models, that is, to adjust a model of

probability for a group of data and to synthesize the result through a distribution

of probability for the parameters of models. The nonlinear models considered to

evaluate the mineralization of organic nitrogen and to illustrate the bayesian

procedure they were: model of Stanford & Smith, model of Marion and model of

Cabrera. The comparison of the models was promoted through the Bayes Factor

(FB) and Bayes Information Criterion (BIC). In this work we had used Gibbs

Sampling and Metropolis Hastings to accomplish inference of the parameters. A

Gibbs Sampling algorithm was implemented on R to get the posterior

distributions of the models parameters. The convergence of the chains was

monitored through graphic analysis, and for the criteria of Geweke and Raftery

& Lewis, implemented in the BOA package, executable in the freeware R. The

model that provided better adjustment quality to the group of data was the model

of Cabrera, being followed by the model of Stanford & Smith and last the one of

Marion. Because of the presented results, it can be attested that the Bayesian

methodology presented good results in the estimate of the parameters of the

models, in other words, the adjustment of models through distributions complete

conditional posteriori constitutes a reliable methodology and accuracy.

Key Words: bayesian inference, nonlinear models, Gibbs Sampling, Metropolis

Hastings, Bayes Factor, Bayesian Information Criterion.

Guidance Committee: Prof. DSc. Joel Augusto Muniz – UFLA (Adviser) and

Prof. DSc. Thelma Sáfadi - UFLA.

1 INTRODUÇÃO

O N potencialmente mineralizável do solo é uma medida da quantidade

de N mineralizado sob condições ótimas de umidade e temperatura, podendo ser

empregado para estimar o nitrogênio (N) que será disponibilizado durante o

ciclo de uma cultura. Durante o período de crescimento das plantas, existe

variação nas quantidades de N inorgânico, principalmente

(nitrato) e

também na capacidade que o solo tem de fornecer N às plantas, por meio da

mineralização. A contribuição relativa dos resíduos e formas orgânicas varia de

acordo com tipo de solo e o sistema de manejo (Stanford & Smith, 1976). Entre

os fatores que contribuem para o incremento da produtividade das culturas, a

disponibilidade de N é um dos mais importantes, pois exerce efeito mais

pronunciado na produção. Apesar disso, o manejo da adubação nitrogenada é

difícil, por ser o N um elemento que apresenta dinâmica complexa e em virtude

do fato da adubação nitrogenada não apresentar efeito residual (Raij, 1991).

N-NO

Compreender as transformações do N no solo constitui um fator

importante para a maximização de seu uso pelas plantas e também para uma

maior eficiência na adubação nitrogenada. Isso porque o nitrogênio é essencial

ao crescimento e à produção das culturas, constituindo, na maioria dos casos, o

elemento que mais limita o crescimento e desenvolvimento das plantas (Silva et

al., 1994).

Na agricultura comercial moderna, grande quantidade de nitrogênio é

adicionada ao solo, constituindo-se fonte potencial de poluição ambiental. A

contaminação dos mananciais com resíduos de

(nitrato) e

(amônia) tem causado sérias preocupações, uma vez que, em excesso (acima de

10 e 0,02 mgL

N-NO

N-NH

-1

), essas formas de nitrogênio podem causar danos à saúde dos

homens e dos animais. Embora tal contaminação seja proveniente de várias

fontes, como resíduos municipais e industriais, as práticas agrícolas, em especial

a adubação do solo com nitrogênio, que têm aumentado abusivamente nos

últimos anos, têm sido consideradas as principais (Ferguson et al., 1991). Em

algumas lavouras, doses elevadas de fertilizantes nitrogenados, quase sempre,

excedem os requerimentos culturais, acarretando eutrofização – lixiviação das

formas de nitrogênio para o lençol freático e ou escoamento superficial das

mesmas em direção aos mananciais (Sainju et al., 1999).

Para que este problema decresça, é necessário que os agricultores usem

adubos nitrogenados e orgânicos em quantidades adequadas que, ao mesmo

tempo, aumentem a produtividade das culturas e diminuam, ao mínimo possível,

os teores de

e no solo (Unlu et al., 1999). De preferência, é

preciso sincronizar melhor a liberação de N pelos resíduos com as fases de maior

requerimento nutricional das culturas. A predição com maior precisão (acurácia)

do N mineralizado possibilita a maior eficiência de uso do nutriente.

N-NO

N-NH

Assim, propõe-se, neste trabalho, desenvolver uma abordagem

bayesiana para predizer as quantidades de N mineralizados por intermédio dos

modelos não lineares, ou seja, ajustar um modelo de probabilidade para um

grupo de dados e sintetizar o resultado por meio de uma distribuição de

probabilidade para os parâmetros dos modelos. Utilizaram-se dados

experimentais (dados reais) associados à mineralização de N em Latossolo do

Sul de Minas Gerais, sob efeito da calagem, descritos em Silva et al. (1994)

(seção 3.1). Na seção (3.5) serão apresentadas as posterioris conjuntas dos

modelos de Cabrera (1993), Marion et al. (1981) e Stanford & Smith (1972). A

comparação dos modelos foi feita por meio do Fator de Bayes e do Critério de

Informação Bayesiano (BIC).

2 REFERENCIAL TEÓRICO

2.1 Comportamento do fenômeno: mineralização do nitrogênio

Segundo Stanford & Smith (1972), a mineralização de N em solo pode

ser descrita por uma equação exponencial simples que relaciona o N

mineralizado, oriundo de um único compartimento (“pool”) de N potencialmente

mineralizável, com o tempo de incubação. Gupta & Reuszer (1967) e Chew et

al. (1976) afirmam que a mineralização da matéria orgânica é maior nos

primeiros períodos de incubação, devido à presença de compostos orgânicos de

fácil decomposição. Black (1968) afirma que fatores, como o manuseio do solo,

a secagem, o congelamento e a moagem, podem acelerar a decomposição da

matéria orgânica.

Com base nessas afirmações e de acordo com Stanford & Smith (1972),

que avaliaram o nitrogênio potencialmente mineralizável do solo (

) por meio

de incubação aeróbica, a mineralização potencial do N orgânico do solo (

)

tem o comportamento conforme descrito na Figura 1.

FIGURA 1: N potencialmente mineralizável (

) com base nas quantidades

acumuladas de N mineralizado durante a incubação sucessiva.

Nota-se, no comportamento dos dados mostrados no gráfico, que, nos

primeiros períodos de incubação, a mineralização do nitrogênio é maior e, à

medida que o tempo decorre, a mineralização da matéria orgânica vai

diminuindo, ou seja, a curva assintótica tende-se a estabilizar. Isso ocorre em

razão de prevalecerem no solo compostos orgânicos que apresentam taxas de

decomposição diferenciadas. De fato, são encontrados no solo compostos

orgânicos nitrogenados facilmente decomponíveis como proteínas, aminoácidos,

açúcares simples, entre outros, e compostos mais resistentes à decomposição,

como é o caso do N associado às substâncias húmicos, açúcares complexos, etc.

Assim, as maiores taxas de decomposição observadas no início da incubação se

explicam em função dos ataques microbianos às moléculas de maior

biodisponibilidade.

À medida que avança o tempo de incubação, diminui a concentração no

solo de moléculas facilmente decomponíveis e aumenta a de compostos

recalcitantes, menos suscetíveis à decomposição. Isso explica o caráter não

linear da mineralização do nitrogênio em função do tempo de incubação, o que

justifica a adoção de modelos não lineares para o estudo do fenômeno

mineralização do nitrogênio.

2.2 Modelos não lineares

Os modelos são classificados como: lineares, não lineares e

linearizáveis. Nos modelos não lineares, não é possível encontrar uma forma

analítica para a estimação dos parâmetros, isto é, as expressões dos estimadores

não apresentam uma solução explícita, exigindo o uso de métodos numéricos

iterativos. A vantagem de se trabalhar com modelos não lineares é que seus

parâmetros possuem uma interpretação biológica e, em muitas situações,

necessitam-se de menos parâmetros nos modelos não lineares, o que simplifica e

facilita a interpretação.

Draper & Smith (1998) classificam os modelos como:

a) modelos lineares: aqueles que são lineares em relação aos parâmetros,

ou seja,

()(

)

XhX

∂

para

e 1, 2, ,ip= … 1, 2, , ;jn= …

em que

p é o número de parâmetros do modelo e n é o número de observações.

Como ilustração, é apresentado o seguinte modelo:

01 2

YXX

ββ ε

++ +

em que

é o erro de caráter aditivo e

, e

ββ

são os parâmetros a serem

estimados. O cálculo das derivadas parciais:

01 2

1, e

YY Y

ββ β

∂∂ ∂

== =

∂∂ ∂

mostra que nenhuma delas depende de algum parâmetro do modelo, portanto, o

modelo é dito linear;

modelos linearizáveis: são modelos não lineares que, por meio de

alguma transformação, se tornam lineares. Seja o modelo,

em que

é o parâmetro a ser estimado e

é o erro de caráter multiplicativo.

Utilizando-se a transformação logarítmica, ou seja, aplicando-se o logaritmo em

ambos os lados da equação, têm-se:

(

)

(

)

() () ()

log log

log log log .

βε

Fazendo

()

(

)

(

)

log ; log ; loggYc e

== =, a equação pode ser escrita como

cX e=+, sendo linear, pois

()

∂

que independe do parâmetro c, mostrando que o modelo original é linearizável.

modelos não lineares: são modelos em que pelo menos uma das

derivadas parciais depende de algum parâmetro do modelo. Seja o

modelo,

βε

em que

são os parâmetros a serem estimados. O cálculo das derivadas

parciais de Y:

ββ

−

∂

∂∂

mostra que a segunda delas depende do parâmetro

, indicando que o modelo

em questão é não linear.

2.2.1 Modelos não lineares adotados

De acordo com Camargo et al. (2002) e Pereira et al. (2005), existem

oito modelos não lineares utilizados para descrever e predizer as relações na

dinâmica da mineralização do nitrogênio no solo. Os oito modelos não lineares

são classificados como: exponencial simples (Cabrera, 1993; Jones, 1984;

Marion et al., 1981 e Stanford & Smith, 1972), exponencial duplo (Inobushi et

al., 1985; Molina et al., 1980), parabólico (Broadbent, 1986) e hiperbólico (Juma

et al., 1984). De forma geral, os modelos exponenciais simples consideram

somente uma forma de N potencialmente mineralizável, decomposto a uma taxa

proporcional à sua concentração. Enquanto os modelos exponenciais duplos

consideram a existência de dois “pools”, ou seja, existem duas taxas ocorrendo

no processo, sendo uma mais estável, característica de compostos mais

resistente, concomitantemente a outra menos estável. A mais estável relacionada

com o N é dificilmente mineralizável e, a menos estável ao N é facilmente

mineralizável (Pereira, 2004).

Três dos modelos estatísticos utilizados são definidos por:

(i) modelo exponencial simples de Stanford & Smith (1972):

(

)

1exp

iii

Nkt=⎡− −⎤+

⎣

⎦

(2.1)

(ii) modelo exponencial simples de Marion et al. (1981):

(

)

1exp

iii

Nkt

⎡

=−−+

⎣

⎤

⎦

(2.2)

(iii) modelo exponencial simples de Cabrera (1993):

(

)

110

1 exp

iiii

Nktkt=⎡− − ⎤+ +

⎣⎦

(2.3)

em que,

é o nitrogênio mineralizado até o tempo

;

é o tempo de incubação, com i = 1, ..., n;

é o nitrogênio potencialmente mineralizável;

é o nitrogênio facilmente mineralizável;

, e kk k

são taxas ou constantes de mineralização;

é um parâmetro que não tem interpretação biológica (constante qualquer);

é um vetor de erros com distribuição normal com média zero e variância

Marion et al. (1981) propuseram um modelo não linear similar ao

modelo descrito por Stanford & Smith (1972) introduzindo um expoente

(

tempo. Cabrera (1993) propôs um modelo não linear que descreve a existência

de dois “pools”, ou seja, duas taxas envolvidas no processo, sendo

a taxa ou

)

constante de mineralização do “pool” mais estável e

a taxa do “pool” menos

estável.

2.3 Inferência bayesiana

A inferência bayesiana nas últimas décadas tem ganhado muitos adeptos

devido ao avanço computacional, ou seja, às técnicas de simulação relativamente

simples, mas extremamente poderosas, que puderam ser implementadas graças

ao avanço dos recursos computacionais. Para Gelman et al. (2003), a inferência

bayesiana é o processo de ajustar um modelo de probabilidade para um grupo de

dados e resumir o resultado por uma distribuição de probabilidade nos

parâmetros do modelo e em quantidades não observáveis como predição para

novas observações.

A estatística é fundamentada em informações a respeito de uma

determinada quantidade de interesse θ e a toda informação relaciona-se um grau

p de incerteza referente a . Um dos objetivos da inferência bayesiana é

minimizar o grau de incerteza, ou seja, reduzir ao mínimo o grau de incerteza a

respeito de (quantidade de interesse). Porém, a incerteza a respeito de θ pode

assumir diferentes graus, que podem ser representados por meio de modelos

probabilísticos para θ . Conseqüentemente, é natural que diferentes

pesquisadores apresentem modelos distintos, pois os mesmos podem apresentar

diferentes graus de “incerteza” sobre θ (Ehlers, 2003).

O grau de “incerteza”, ou a informação que se tem a respeito de , é

representado probabilisticamente por meio de

( )

. A informação sobre θ pode

ser aumentada observando-se uma quantidade aleatória Y relacionada com θ , na

qual a distribuição amostral define esta relação. Uma vez observada

, a quantidade de informação sobre θ aumenta e o teorema de Bayes é a

( | )

Yy=

regra de atualização utilizada para mensurar este aumento de informação

(Ehlers, 2003),

()

(

)

()

(

)

(

)

()

(

)

(

)

()

,| |

py pyp py p

py py

θθθθ

θθ

== =

∫

. (2.4)

Observe que, em (2.4),

()

1 py

desempenha o papel de uma constante

normalizadora de

(

)

|pyθ

, pois não depende de . Para um valor fixo de y, a

função

(

)

|pyθ

fornece a verossimilhança de cada um dos possíveis valores de

, ou seja,

, enquanto

() (

|py L yθθ=

)

;

( )

é chamada de distribuição

priori de . Essas duas fontes de informação, priori e verossimilhança, são

combinadas, levando à distribuição a posteriori de , . Portanto, a forma

usual do teorema de Bayes é:

θ ( | )

()()(

)

yLypθθ∝ x θ (2.5)

É importante ressaltar que essa forma simplificada do teorema de Bayes

dada em (2.5) é útil para a resolução de problemas que envolvam estimação de

parâmetros, já que o denominador é apenas uma constante normalizadora.

Porém, em outras situações, como seleção de modelos, este termo tem um papel

crucial (Ehlers, 2003).

As inferências sobre a distribuição a posteriori conjunta consistem na

verdade em encontrar uma distribuição para um parâmetro específico,

. A

forma de encontrar tal distribuição, chamada distribuição marginal de

, é

integrando a distribuição posteriori conjunta em relação aos outros parâmetros

do modelo, ou seja,

(

)

(

)

iiii

θθ θ

−

∫∫



(2.6)

em que

é o conjunto complementar de parâmetros

para

(

)

111

,, , ,,

iii

θθθθ θ

−−+

= ……

(Rosa, 1998).

Porém, muita das vezes, a forma analítica da marginal (2.6) é complexa

e multidimensional ou, até mesmo, impossível de ser calculada analiticamente,

ou seja, a inferência exata somente será possível se estas integrais puderem ser

calculadas analiticamente, caso contrário devem-se usar aproximações. Nesse

caso, métodos numéricos iterativos, como método de Monte Carlo via cadeias de

Markov (MCMC) são utilizados para gerar valores de uma distribuição

condicional a posteriori para cada parâmetro (Sorensen, 1996). Dentre os

métodos de simulação que fazem o uso de cadeias de Markov destacam-se o

amostrador de Gibbs e o Metropolis Hastings (Hastings, 1970).

O amostrador de Gibbs é uma alternativa para a obtenção de (2.6),

método inicialmente utilizado por Gelfand & Smith (1990), na estatística, o qual

foi elaborado por Geman & Geman (1984). Gamerman (1997) também utiliza o

método MCMC para amostragem da densidade a posteriori e, para os casos em

que não existem formas específicas para as distribuições condicionais completas,

a amostragem Gibbs pode ser implementada com o auxílio do algoritmo

Metropolis & Hastings (Hastings, 1970; Metropolis et al., 1953), para realizar o

sorteio na condicional completa da amostragem Gibbs.

A motivação para o uso da amostragem Gibbs é que a real distribuição

pode ser aproximada por uma distribuição empírica de n valores, de tal forma

que n seja grande o suficiente para que a amostragem Gibbs atinja a

convergência (Zeger & Karin, 1991). Gelfand & Smith (1990) utilizam a mesma

idéia ilustrando a inferência bayesiana em um conjunto de dados normais,

usando o amostrador de Gibbs. Casella & George (1992) utilizam o amostrador

de Gibbs para ilustrar o fato de que as condicionais completas, quando

determinam as marginais, isto é, quando se tem uma amostra suficientemente

n →∞

grande de um determinado parâmetro, dados os demais, obtém-se uma

distribuição empírica que se aproxima suficientemente da distribuição marginal.

2.3.1 Obtenção das distribuições posterioris condicionais completas

Em geral, cada uma das componentes

pode ser uni ou

multidimensional, a distribuição condicional completa é a distribuição da i-

ésima componente de condicionada em todas as outras componentes. A

condicional completa pode ser obtida a partir da distribuição conjunta da

seguinte forma:

()

θθ

−

∫

(2.7)

em que

, e

(

)

111

,, , ,,

iiij

θθθθ θ

−−+

= ……

(

)

,,,,

θθ θ θ= ……

()

θ é a

distribuição conjunta, sendo j o número de parâmetros do modelo. Assim, para

se obter a distribuição condicional completa de

, basta pegar os termos da

distribuição conjunta que não dependem de

(Sorensen, 1996 e Ehlers, 2003).

2.3.2 Priori não informativa

A priori é caracterizada como não informativa quando a verossimilhança

é dominante ou quando se deseja representar o desconhecimento sobre θ . A

primeira idéia é pensar em todos os possíveis valores para como igualmente

prováveis, isto é,

()

vθ ∝ , com v sendo uma constante e limitado em um

intervalo real. Muitas vezes, a priori, é a forma de se quantificar a incerteza

sobre o parâmetro, porém, a priori não-informativa não representa

necessariamente o desconhecimento do pesquisador sobre o parâmetro, mas

deve ser usada também de forma a viabilizar a inferência a posteriori (Box &

Tiao, 1992).

2.3.3 Intervalo de credibilidade

Na inferência bayesiana, toda informação probabilística que se tem a

respeito de um determinado parâmetro é representada por meio de sua

distribuição a posteriori. A partir de tal distribuição, estimam-se parâmetros

como média, mediana e moda a posteriori, baseando-se nos valores observados

da distribuição a posteriori de .

No entanto, essa informação restringe-se a um único valor, ou seja, toda

informação presente na distribuição a posteriori de é condensada em uma

estimação pontual. Portanto, fica evidente a necessidade de associarmos uma

medida da precisão de tal estimativa (média, mediana ou moda). Para isso, é

necessário incorporar o conceito de intervalo de credibilidade que permite

mensurar a precisão com que este número foi estimado.

Define-se

[

]

,Cab= como intervalo de credibilidade de ,

ou nível de credibilidade

(

, para se

()

100 1 %α−

)

1 α− θ

[

]

1PCθ∈≥−α. Note que a

definição expressa, de forma probabilística, a pertinência ou não de θ ao

intervalo. Assim, quanto menor for o tamanho do intervalo, mais concentrada é a

distribuição do parâmetro, ou seja, o tamanho do intervalo informa sobre a

dispersão do parâmetro (Ehlers, 2003). θ

2.3.3.1 Intervalo de máxima densidade a posteriori ou intervalo HPD

A sigla HPD corresponde à abreviação da expressão “Highest Posterior

Density”, ou seja, região de alta ou máxima densidade a posteriori. A estimação

de intervalo de densidade de probabilidade máxima (HPD) é um meio de gerar

intervalos posterioris bayesianos. Paulino et al. (2003) afirmam que existe uma

infinidade de regiões de credibilidade com o mesmo grau de credibilidade.

Portanto, o intervalo HPD é um intervalo de credibilidade de comprimento

mínimo, obtido a partir dos valores de com maior densidade a posteriori. No

pacote BOA, executável no “freeware” R, utiliza-se o algoritmo descrito por

Chen & Shao (1999) para estimar os intervalos HPD, sob a suposição de

distribuições posterioris marginais unimodal.

2.3.4 Erro de Monte Carlo

No contexto de inferência bayesiana, o método de Monte Carlo é

utilizado para o cálculo de integrais definidas, utilizando-se aproximações.

De acordo com Ehlers (2003), o estimador de Monte Carlo é dado por:

()

∑

em que é uma aproximação para a integral que se deseja calcular, portanto,

precisamos estudar o erro . Uma vez que as gerações são independentes,

então, pela Lei Forte dos Grandes Números, segue que converge quase

certamente para , isto é,

II−

() ()

XEgX

→∞

⎡

⎤

→

⎣

⎦

∑

Além disso, definindo

()

Var

Xσ

⎡

⎤

⎣

⎦

e assumindo que esta variância

existe o erro padrão de Monte Carlo é uma estimativa consistente de σ , dada

por:

()

⎡

⎤

= −

⎣

⎦

∑

isto é, a aproximação pode ser tão acurada quanto se deseje, bastando aumentar

o valor de n (número de iterações).

Este método consiste em gerar vários números aleatórios uniformemente

distribuídos no intervalo de interesse, avaliando-se o valor da função a ser

integrada em cada um destes pontos e, finalmente, calculando-se a média desses

valores obtidos para se obter a aproximação da integral definida. Como o

método de Monte Carlo é uma aproximação, não se pode afirmar que o valor

encontrado por ele é o valor exato da integral. Entretanto, quanto mais pontos

forem gerados, maior será a confiabilidade na aproximação obtida.

3 MATERIAL E MÉTODOS

3.1 Material

Para ilustração do estudo dos três modelos, utilizaram-se dados

experimentais associados à mineralização de N em Latossolo do Sul de Minas

Gerais sob efeito da calagem, descritos em Pereira et al. (2005) e Silva et al.

(1994). Assim, foram consideradas as quantidades acumuladas de N

mineralizado durante onze tempos de incubação (1, 2, 3, 4, 6, 8, 10, 12, 14, 21 e

28 semanas). A acidez do solo foi corrigida visando-se elevar a saturação por

bases a 60%, com o pH em água do solo atingindo valores na faixa de 6,0-6,2. A

extração do N mineral (nitrato + amônio) formado durante a incubação foi feita

pela lavagem periódica do solo com solução de CaCl

(Cloreto de Cálcio) 0,01

mol L

-1

. A quantificação do N mineralizado foi feita em destilador de arraste de

vapores, sendo o N-amônio quantificado após a adição de MgO (Óxido de

Magnésio) no extrato de CaCl

e o nitrato, após o uso de liga de Devarda no

extrato remanescente. As amostras de solo foram incubadas em condições

ambientais do Laboratório de Relação Solo-Planta do DCS/UFLA

(Departamento de Ciências do Solo da Universidade Federal de Lavras), com a

temperatura ambiente oscilando entre 21

C e 28

C. O delineamento

experimental utilizado foi o inteiramente ao acaso, com três repetições.

3.2 Métodos

Para realizar a inferência sobre os parâmetros é necessário determinar as

marginais a posteriori. De posse das marginais, utiliza-se o método MCMC para

realizar amostragem das densidades a posterioris e, para os casos em que não

existem formas específicas para as distribuições condicionais completas, a

amostragem Gibbs é implementada por intermédio do algoritmo Metropolis

Hastings. Serão apresentados alguns métodos de reamostragem com o intuito de

facilitar a compreensão do amostrador de Gibbs e a amostragem por rejeição

com algoritmo Metropolis Hastings.

3.2.1 Métodos de reamostragem

A idéia dos métodos de reamostragem é gerar valores em duas etapas.

Na primeira etapa, geram-se valores de uma distribuição auxiliar conhecida no

sentido de aproximar-se da distribuição de interesse, ou seja, da posteriori. Na

segunda etapa, utiliza-se um mecanismo de correção a fim de que os valores

amostrados da distribuição auxiliar sejam representativos (ao menos

aproximadamente) da distribuição a posteriori. O mecanismo de correção não

necessariamente é estocástico. A densidade de interesse será denotada por

(

)

xθ e a densidade auxiliar é denotada por (Gamerman, 1997).

()

q θ

3.2.1.1 Método da rejeição

Considere uma densidade auxiliar , da qual sabemos gerar

amostras. A idéia é utilizar

para fazer geração de

()

q θ

()

q θ

(

)

xθ . A única

restrição matemática sobre

é que exista uma constante , tal que

, para todo . O método gera valores da densidade

proporcional a

()

q θ

A< ∞

() (

|pxAqθ<

)

(

xθ , se apenas o núcleo de

(

)

xθ é conhecido. Basta que

sirva como envelope, cobrindo todo o gráfico de

()

q θ

(

)

xθ , independente

da escala utilizada (Gamerman, 1997).

O método consiste em gerar um valor

de uma distribuição auxiliar

e aceitar este valor como sendo da distribuição a posteriori com

probabilidade

()

q θ

() ()

xAqθθ

)

. Caso contrário, não é aceito como um valor

gerado da posteriori e o processo é repetido até que um valor seja aceito. O

método também funciona se em vez da posteriori,

(

xθ , que, em geral, é

desconhecida, utilizarmos a sua versão não normalizada, ou seja,

(Ehlers, 2003).

()(

|px pθθ

)

Em algumas situações, toma-se a priori

(

)

θ como sendo a densidade

auxiliar, ou seja,

, então, a constante A deve ser tal que

() ()

qpθ= θ

()

xθ < A

)

. Esta desigualdade só pode ser satisfeita se tomarmos A como

sendo o valor máximo da função de verossimilhança, isto é,

, em

que

é o estimador de máxima verossimilhança de θ . Então, a probabilidade

de aceitação de

se simplifica em

(

|Apxθ=



()

|px pxθ



)

|θ , que é sempre menor ou

igual a 1.

Pode-se utilizar o seguinte algoritmo para gerar valores da posteriori:

1) gerar um valor

da distribuição auxiliar (conhecida);

()

q θ

2) gerar

;

()

0,1uU∼

3) aceitar

como valor da posteriori se

() (

|up xAqθ< θ, caso contrário

rejeitar

e retornar ao item 1.

Um problema técnico associado ao método é a necessidade de se

maximizar a função de verossimilhança, o que pode não ser uma tarefa simples

em modelos mais complexos. Se este for o caso, então, o método de rejeição

perde o seu principal atrativo, que é a simplicidade (Gamerman, 1997).

3.2.1.2 Método de reamostragem ponderada

O método de reamostragem ponderada tem o mesmo princípio do

método da rejeição, ou seja, parte de valores gerados de uma distribuição

auxiliar

fácil de ser amostrada, mas com a vantagem de que ela não

precisa ser um envelope para

()

q θ

(

)

xθ , isto é, não há necessidade de

maximização da verossimilhança. A desvantagem do método é a de que ele

fornece valores apenas aproximadamente distribuídos segundo a posteriori

(Gamerman, 1997).

Suponha que temos uma amostra

gerada de uma distribuição

auxiliar

e, a partir dela, construímos pesos com , em que

θ… θ

()

q θ

1, ,i= …

()()

pxq

θθ

∑

O método consiste em tomar uma segunda amostra (reamostra) de

tamanho m da distribuição discreta em

, com probabilidades

. Aqui, também, não é necessário que se conheça completamente a

posteriori, mas apenas o produto da priori pela verossimilhança (versão não

normalizada), já que, neste caso, os pesos não se alteram (Ehlers, 2003).

θ… θ

ww…

Tomando novamente a priori como densidade auxiliar, isto é,

, os pesos se simplificam para:

() ()

qpθ=

()

∑

, .

1, ,in= …

O algoritmo para gerar valores (aproximadamente) da posteriori é:

1) gerar valores

de uma distribuição auxiliar (conhecida);

θθ…

()

q θ

2) calcular os pesos

, .

1, ,in= …

3) reamostrar valores com probabilidades

ww…

3.2.2 Método de Monte Carlo via cadeias de Markov (MCMC)

Todos os métodos de simulação apresentados até agora geram os valores

de forma independente (não-Markovianos), isto é, as observações sucessivas são

estatisticamente independentes. O método de Monte Carlo via cadeias de

Markov (MCMC) consiste de simulações iterativas, baseadas em cadeias de

Markov e, conseqüentemente, os valores gerados não serão mais independentes

(Gamerman, 1997). Os métodos MCMC mais utilizados são o amostrador de

Gibbs e o algoritmo de Metropolis-Hastings.

3.2.2.1 Amostrador de Gibbs

No amostrador de Gibbs, a cadeia irá sempre se mover para um novo

valor, isto é, não existe mecanismo de aceitação-rejeição. As transições de um

estado para outro são feitas de acordo com as distribuições condicionais

completas

(

)

p θθ

−

, em que , sendo j o número de

parâmetros do modelo (Ehlers, 2003).

(

)

111

,, , ,,

iii

θθθθ θ

−−+

= ……

O amostrador de Gibbs é um processo Markoviano, por sustentar uma

memória curta, sendo a variável gerada dependente apenas da variável gerada

anteriormente. O algoritmo de Gibbs fornece uma forma alternativa para gerar

sucessivos valores de uma distribuição condicional completa. Segundo

Gamerman (1997), o amostrador de Gibbs pode ser descrito da seguinte forma:

1) inicie o contador de iterações da cadeia l = 1 e escolha valores iniciais

;

() ()

()

θθ θ= …

2) obtenha um novo valor

, a partir de

() () ()

(

ll l

θθ θ= …

)

(

)

−

, por meio de

sucessivas gerações de valores

(

)

(

)

(

)

(

)

(

)

() () ( ) ( )

()

() () () ()

(

)

11 1

1123

2213

12 1

|,,,

|, ,,

|,,,

lll

jj j

θθθθθ

−− −

−−

−

∼…



∼…

3) mude o contador l para l +1 e retorne ao passo 2, até atingir a convergência.

3.2.2.2 Algoritmo de Metropolis-Hastings

Em alguns casos, a forma da distribuição condicional completa não é

reconhecida, o que inviabiliza sua amostragem pelos métodos convencionais.

Alguns métodos adicionais podem ser utilizados, como, por exemplo,

amostragem por rejeição com o algoritmo de Metropolis-Hastings.

O algoritmo de Metropolis-Hastings utiliza a mesma idéia dos métodos

de rejeição vistos anteriormente, ou seja, um valor é gerado da distribuição

auxiliar e aceito com uma dada probabilidade. Esse mecanismo de correção

garante a convergência da cadeia para a distribuição de equilíbrio (Ehlers, 2003).

Suponha que a cadeia esteja no estado θ e um valor

é gerado de uma

distribuição auxiliar proposta

. Note que a distribuição auxiliar proposta

pode depender do estado atual da cadeia. O novo valor

é aceito com a

seguinte probabilidade:

()

.|q θ

()

()()

()

,min1,

pxq

θθθ

αθθ

θθθ

⎛⎞

⎟

⎜

⎟

⎜

⎟

⎜

⎟

⎜

⎟

⎜

⎝⎠

(3.1)

em que p é a distribuição de interesse.

É notório que a cadeia pode permanecer no mesmo estado por muitas

iterações e, na prática, costuma-se monitorar isso calculando-se a porcentagem

média de iterações para as quais novos valores são aceitos. O algoritmo de

Metropolis-Hastings pode ser descrito da seguinte forma:

1) inicie o contador de iterações da cadeia t = 0 e escolha valores iniciais

;

2) gere um novo valor

da distribuição ;

()

.|q θ

3) calcule a probabilidade de aceitação

(

)

,αθθ por meio de (3.1) e gere

;

()

0,1uU∼

4) se

então aceita-se o novo valor e faça , caso contrário,

rejeite e faça

;

α≤

()

= θ

() ()

1tt

θθ

5) incremente o contador de t para t + 1 e volte ao passo 2.

Maiores detalhes sobre este algoritmo podem ser encontradas em Chib

& Greenberg (1995), Gamerman (1997) e Hastings (1970).

3.3 Diagnósticos de convergência

Uma questão de ordem prática é como os valores iniciais influenciam o

comportamento da cadeia. Considere uma longa realização da cadeia com

comprimento (número de iterações) t = l + gm , em que l é o número inicial de

iterações necessárias para que a cadeia atinja o estado de equilíbrio, m é o

número das iterações que serão utilizadas na aplicação do método de Monte

Carlo e g o espaçamento (pulo) entre iterações sucessivas, delineado para

eliminar a autocorrelação entre as consecutivas iterações. O período constituído

pelas l primeiras iterações costuma designar-se por período de aquecimento

(burn-in). Este período pode ser mais ou menos longo e, para contornar este

“problema”, há métodos de diagnóstico para verificar se a cadeia se encontra ou

não no estado de equilíbrio (Paulino et al., 2003).

Conforme o número de iterações aumenta, a cadeia gradualmente

esquece os valores iniciais e, eventualmente, converge para uma distribuição de

equilíbrio. Assim, em aplicações práticas, é comum que as iterações iniciais

sejam descartadas, como se formassem uma amostra de aquecimento.

Os métodos ou critérios utilizados para diagnosticar, ou seja, para

monitorar a convergência, podem ser classificados como informais e formais. Os

métodos informais foram inicialmente propostos por Gelfand & Smith (1990).

Esses autores sugeriram técnicas gráficas para a verificação da convergência,

como, por exemplo, a técnica gráfica que consiste em observar a trajetória da

cadeia ao longo das iterações. Se o gráfico, após um período inicial, apresenta

repetidamente o mesmo comportamento qualitativo e quantitativo, então, pode-

se concluir pela convergência da cadeia (Gamerman, 1997).

Dentre os métodos formais, podem-se citar os seguintes critérios de

convergência: Gelman & Rubin (1992), Geweke (1992), Heidelberger & Welch

(1983), Raftery & Lewis (1992) e outros mais. Esses métodos e outros foram

comparados no trabalho de Cowles & Carlin (1996), que concluíram que não se

pode afirmar qual deles é o mais eficiente.

No presente trabalho, a monitoração da convergência das cadeias foi

realizada por meio de análise gráfica, critério de Geweke (1992) e critério de

Raftery & Lewis (1992). Os critérios de Gelman & Rubin, Raftery & Lewis,

Geweke, Heidelberger & Welch e outros estão implementados no pacote BOA,

executável no freeware R, conforme R Development Core Team (2004).

3.3.1 Critério de Raftery & Lewis

O critério proposto por Raftery & Lewis (1992) é um método que estima

quantas iterações são necessárias para que o amostrador de Gibbs apresente

convergência à distribuição estacionária, propondo quantidades necessárias para

isso. O método fornece as estimativas do “burn-in”, que é o número de iterações

que devem ser descartadas, o número de iterações que devem ser computadas N

(Total) e a distância mínima de uma iteração à outra, para se obter a subamostra

aproximadamente independente (“thin”). Ou seja, o método fornece o tamanho

do descarte que deve ser dado para extrair o efeito dos valores iniciais e o salto

de uma iteração para a outra, para retirar o efeito de dependência entre os valores

amostrados. A regra de decisão do critério de Raftery & Lewis é baseada no

fator de dependência. O fator de dependência é responsável pelo acréscimo

multiplicativo ao número de iterações necessárias para se alcançar a

convergência, devido à autocorrelação dentro da seqüência. De acordo com

Raftery & Lewis (1992), se o fator de dependência for maior que 5, pode-se

concluir que a cadeia não atingiu convergência.

3.3.2 - Critério de Geweke

O critério proposto por Geweke (1992) consiste em indicar a

convergência da média a posteriori da quantidade amostrada,

(

)

, em uma

única cadeia, isto é, os valores

(

)

(

)

são calculados a cada iteração da

simulação, formando uma série temporal. A partir da cadeia, a variância

assintótica

da medida de

()

(

)

pode ser estimada.

O critério de Geweke consiste nos seguintes passos:

Passo 1: divide-se a cadeia de tamanho N em duas seqüências. A primeira

seqüência possui primeiras iterações e a segunda restantes iterações.

Geweke (1992) sugere que as médias devem ser estimadas após o descarte de

algumas iterações iniciais e que sejam considerados

0,1 e 0,5

NNNN

Passo 2: estimam-se as médias

ˆˆ

e as variâncias assintóticas

(

)

(

)

ˆˆ

0 e 0

, respectivamente. As variâncias são determinadas pela

estimação da densidade espectral.

Pode-se mostrar que, se as razões

NN NN

são fixas e ,

então:

N →∞

() ()

()

ˆˆ

0,1

ˆˆ

−

∼

Portanto, um teste pode ser construído e se a diferença padronizada entre

as médias for muito grande, existe indicação de ausência de convergência. A

regra de decisão do critério de Geweke é baseada na estimativa do p-valor. De

acordo com Geweke (1992), se o p-valor for menor que um nível de

significância adotado pelo pesquisador, conclui-se que a cadeia não atingiu

convergência.

3.4 Comparação dos modelos não lineares

A comparação dos modelos foi realizada por intermédio do Fator de

Bayes e do Critério de Informação Bayesiano (BIC).

3.4.1 Fator de Bayes (FB)

Suponha dois modelos

. Segundo Berkhof et al. (2003), Han &

Carlin (2001) e Paulino et al. (2003), o Fator de Bayes está formalmente

definido como a razão entre as posterioris e as prioris, ou seja,

()()

()

pM x pM x

pM pM

⎡

⎤

⎣⎦

⎡⎤

⎣

⎦

em que

(

)

M é a probabilidade a priori e

()()()()

iii

Mx pMpxM px

⎡

⎤

⎣

⎦

é a probabilidade a posteriori do modelo

com . Observe que se as

prioris forem iguais, então

1, 2i =

é maior que 1 se

apresentar uma

probabilidade posteriori maior, ou seja, o modelo 1 é mais plausível ou

verossímil do que o modelo 2, então, opta-se por escolher o modelo 1 para

predizer o fenômeno em questão.

3.4.2 Critério de Informação de Bayesiano (BIC)

De acordo com Carlin & Louis (2000), o Critério de Informação

Bayesiano (BIC) pode ser definido da seguinte forma:

()

{

}

()

2ln |, ln

iiii

IC E L x M p nθ

⎡

⎤

= −

⎣

⎦

em que n é o tamanho da amostra.

A regra de decisão baseada no BIC consiste na escolha do modelo que

apresentar maior valor de BIC

3.5 Posteriori conjunta

3.5.1 Posteriori conjunta para o modelo não linear de Stanford & Smith

Para a obtenção da posteriori conjunta (3.5) do modelo não linear de

Stanford & Smith (2.1), foram estabelecidas:

i) uma priori não informativa (seção 2.3.2) para os parâmetros

, ou seja,

e Nk

()

Nk w∝ . (3.2)

Enquanto que, para a precisão

, , assume-se como priori:

()

τσ

−

ii) uma distribuição Gama com parâmetros

()

, ,

αβ

()

{

}

|, exp

eee e ee

ταβ τ βτ

−

∝−. (3.3)

Supondo que

, então, a verossimilhança para o

modelo (2.1) pode ser escrita da seguinte forma:

(

eiidNτ

−

∼

)

() ()

{}

() ()

{}

() ()

{}

| , , , exp 1 exp

|,,, exp 1exp

eei

ee i i

LyN k H y N kt

πσ

ττ

⎧⎫

⎪⎪

⎡⎤

= −−−−

⎨⎬

⎣⎦

⎪⎪

⎩⎭

⎧⎫

⎪⎪

⎡

⎤

= −−−−

⎨⎬

⎣

⎦

⎪⎪

⎩⎭

⎧⎫

⎪⎪

⎡⎤

∝−−−−

⎨⎬

⎣⎦

⎪⎪

⎩⎭

∏

() ()

{}

| , , , exp 1 exp

ee i i

LyN k H y N kt

ττ

⎧⎫

⎪⎪

⎡

⎤

∝− −−−

⎨⎬

⎣

⎦

⎪⎪

⎩⎭

∑

(3.4)

em que i = 1, ..., n ,

{

}

,,,= …

yyy ye H são os hiperparâmetros: .

αβ

Agora de posse da verossimilhança (3.4), das prioris (3.2) e (3.3) e por

intermédio do teorema de Bayes (2.5) a distribuição posteriori conjunta do

modelo (2.1) pode ser descrita da seguinte forma:

(

)

(

)

(

)

(

)

(

)

(

)

(

)

000

,, |, | ,, , ,, |

,, |, | ,, , , |

pN k yH LyN k H pN k H

Nk yH LyNk HpNkp H

ττ

∝

() ()

{}

{

}

,, |, exp 1 exp

exp

ee i i

eee

pN k yH y N kt

ττ

τβτ

−

⎧⎫

⎪⎪

⎡⎤

∝− −−−

⎨⎬

⎣⎦

⎪⎪

⎩⎭

−

∑

x x

()

{}

-1

,, |,

exp 2 1 exp

ei i

pN k yH

yN kt

ττ

∝

⎧⎫

⎡

⎤

⎪⎪

⎡

⎤

⎢

⎥

− + −−−

⎨⎬

⎣

⎦

⎢

⎥

⎪⎪

⎣

⎦

⎪⎪

⎩⎭

∑

. (3.5)

3.5.2 Posteriori conjunta para o modelo não linear de Marion

De modo análogo, para a obtenção da posteriori conjunta (3.9) do

modelo não linear de Marion (2.2), foi estabelecida:

i) uma priori não informativa para os parâmetros

, e Nkb

()

Nkb v∝ . (3.6)

Enquanto que, para a precisão

, , assume-se como priori:

()

τσ

−

ii) uma distribuição Gama com parâmetros

()

, ,

αβ

()

{

}

|, exp

eee e ee

ταβ τ βτ

−

∝−. (3.7)

Dado que

, então, a verossimilhança para o modelo

(2.2) pode ser escrita da seguinte forma:

(

eiidNτ

−

∼

)

()

{

}

()

{}

()

{}

| , , , , exp 1 exp

eei

ee i i

LyN kb H y N kt

πσ

ττ

⎧⎫

⎪⎪

⎡⎤

= −−−−

⎨⎬

⎢⎥

⎣⎦

⎪⎪

⎩⎭

⎧⎫

⎪⎪

⎡

⎤

= −−−−

⎨⎬

⎢

⎥

⎣

⎦

⎪⎪

⎩⎭

⎧⎫

⎪⎪

⎡⎤

∝−−−−

⎨⎬

⎢⎥

⎣⎦

⎪⎪

⎪

⎩⎭

∏

⎪

∏

()

{

}

| , , , , exp 1 exp .

ee i i

LyN kb H y N kt

ττ

⎧⎫

⎪⎪

⎡

⎤

∝− −−−

⎨⎬

⎢

⎥

⎣

⎦

⎪⎪

⎩⎭

∑

(3.8)

em que i = 1, ..., n,

{

}

,,,= …

yyy y e H são os hiperparâmetros: .

αβ

De posse da verossimilhança (3.8), das prioris (3.6) e (3.7) e por

intermédio do teorema de Bayes (2.5), a distribuição posteriori conjunta do

modelo (2.2) pode ser descrita da seguinte forma:

(

)

(

)

(

)

(

)

(

)

(

)

(

)

000

,,, |, | ,,, , ,,, |

,,, |, | ,,, , ,, |

pN kb yH LyN kb H pN kb H

Nkb yH LyNkb HpNkbp H

ττ

∝

()

{}

{

}

,,, |, exp 1 exp

exp

ee i i

eee

pN kb yH y N kt

ττ

τβτ

−

⎧⎫

⎪⎪

⎡⎤

∝− −−−

⎨⎬

⎢⎥

⎣⎦

⎪⎪

⎩⎭

−

∑

x x

()

{}

,,, |,

exp 2 1 exp .

ei i

pN kb yH

yN kt

ττ

+ −

∝

⎧⎫

⎡

⎤

⎪⎪

⎡⎤

⎢

⎥

− + −−−

⎨⎬

⎢⎥

⎣⎦

⎢

⎥

⎪⎪

⎣

⎦

⎪⎪

⎩⎭

∑

(3.9)

3.5.3 Posteriori conjunta para o modelo não linear de Cabrera

Analogamente, para a obtenção da posteriori conjunta (3.13) do modelo

não linear de Cabrera (2.3), foram estabelecidas:

i) uma priori não informativa para os parâmetros

11 0

, e Nk k

()

110

Nkk z∝ . (3.10)

Para a precisão

, , assume-se como priori:

()

τσ

−

ii) uma distribuição Gama com parâmetros

()

, ,

αβ

()

{

}

|, exp

eee e ee

ταβ τ βτ

−

∝−. (3.11)

Uma vez que

, então, a verossimilhança para o modelo

(2.3) pode ser escrita da seguinte forma:

(

eiidNτ

−

∼

)

() ()

{}

() ()

{}

() ()

110 1 1 0

110 1 1

|,,,, exp 1exp

eei i

ee i i

LyNkk H y N kt kt

LyNkk H y N kt

πσ

ττ

⎧⎫

⎪⎪

⎡⎤

= −−−−−

⎨⎬

⎣⎦

⎪⎪

⎩⎭

⎧⎫

⎪⎪

⎡⎤

= −−−−−

⎨⎬

⎣⎦

⎪⎪

⎩⎭

⎡⎤

∝−−−−

⎣

∏

{}

0 i

⎧⎫

⎪⎪

−

⎨⎬

⎦

⎪⎪

⎩⎭

() ()

{}

110 1 1 0

| , , , , exp 1 exp .

ee i ii

LyNkk H y N kt kt

ττ

⎧⎫

⎪⎪

⎡⎤

∝− −−−−

⎨⎬

⎣⎦

⎪⎪

⎩⎭

∑

(3.12)

em que i = 1, ..., n,

{

}

,,,= …

yyy y e H são os hiperparâmetros: .

αβ

De posse da verossimilhança (3.12), das prioris (3.10) e (3.11) e por

intermédio do teorema de Bayes (2.5), a distribuição posteriori conjunta do

modelo (2.3) pode ser descrita como:

()()(

()()()

110 110 110

,,, |, | ,,,, ,,, |

,,, |, | ,,,, ,, |

pNkk yH LyNkk HpNkk H

)

()

Nkk yH LyNkk HpNkk p H

ττ

∝

() ()

{}

{

}

110 1 1 0

,,, |, exp 1exp

exp

ee i ii

eee

pN kk yH y N kt kt

ττ

τβτ

−

⎧⎫

⎪⎪

⎡⎤

∝− −−−−

⎨⎬

⎣⎦

⎪⎪

⎩⎭

−

∑

x x

()

{}

110

,,, |,

exp 2 1 exp .

ei ii

pNkk yH

yN kt kt

ττ

+ −

∝

⎧⎫

⎡

⎤

⎪⎪

⎡⎤

⎢

⎥

− + −−−−

⎨⎬

⎣⎦

⎢

⎥

⎪⎪

⎣

⎦

⎪⎪

⎩⎭

∑

(3.13)

A estimação dos parâmetros foi feita por meio de rotinas elaboradas no

freeware R e, para o procedimento de amostragem, utilizou-se o amostrador de

Gibbs e o Metropolis Hastings. No processo, foram geradas 450.000 iterações,

considerando-se um “burn-in” de 50.000, isto é, as 50.000 primeiras iterações

foram descartadas, e um “thin” a cada 100 iterações para assegurar a

independência da amostra. Os modelos não lineares (2.1), (2.2) e (2.3) foram

comparados por intermédio do Fator de Bayes e do BIC (Paulino et al., 2003).

Dessa forma, obteve-se uma amostra de tamanho 4.000 para cada cadeia

(parâmetro). O pacote Bayesian Output Analysis (BOA) do Sistema

Computacional Estatístico R, conforme R Development Core Team (2004), foi

utilizado para realizar inferências sobre os parâmetros dos modelos (2.1), (2.2) e

(2.3).

4 RESULTADOS E DISCUSSÃO

4.1 Obtenção das posterioris condicionais completas para o modelo de

Stanford & Smith

De acordo com a equação (2.7), na seção (2.3.1), a posteriori

condicional completa para o parâmetro

do modelo (2.1), segundo a sua

posteriori conjunta (3.5), é dada por:

()

,, |,

|, ,,

,, |,

pN k yH

Nk yHdN

∫

() ()

{}

|, ,, exp 1 exp .

pN k yH y N kt

⎧⎫

⎪⎪

⎡

∝− −−−

⎨⎬

⎣

⎪⎪

⎩⎭

∑

⎤

⎦

(4.1)

Considerando-se

e como vetor dos elementos

para o tempo i, com

1 exp( )= −

1, 2, ,=…in

{

}

,,,= …

yyy y pode-se definir (4.1)

como:

()()()

[]

() () ()

000

00 00

000

|, ,, exp

exp ' ' ' '

exp ' ' ' ' '

pN k yH y Na y Na

yy yNa aN y aN Na

yy yNa aNy Naa

NaaayaaNaaay

−−

⎧⎫

⎪⎪

∝−− −

⎨⎬

⎪⎪

⎩⎭

⎧⎫

⎪⎪

∝− − − +

⎨⎬

⎪⎪

⎩⎭

⎧⎫

⎪⎪

⎡⎤

∝− − − +

⎨⎬

⎢⎥

⎣⎦

⎪⎪

⎩⎭

⎧⎫

⎪⎪

⎡

⎤⎡

∝− − −

⎨⎬

⎤

⎢

⎥⎢ ⎥

⎣

⎦⎣

⎪

⎩⎭

⎪

⎦

que é o núcleo de uma distribuição normal. Portanto, a posteriori condicional

completa para

é uma normal, ou seja,

(

)

(

)

(

)

(

)

|, ,, ' ', '

Nk yH Naa ay aaτ

−

∼

−

, (4.2)

em que

()

'1exp

⎡⎤

= −−

⎣⎦

∑

aa kt

()

1exp

−

⎡

⎤

−−

⎣

⎦

⎡

⎤

−−

⎣

⎦

∑

kt y

aa ay

Logo, a posteriori condicional completa para

pode ser reescrita da seguinte

forma:

()

1exp

|, ,, , .

1exp 1exp

kt y

Nk yH N

kt kt

⎛⎞

⎡⎤

−−

⎟

⎜

⎟

⎣⎦

⎜

⎟

⎜

⎟

⎜

⎟

⎜

⎟

⎜

⎟

⎡

⎤⎡ ⎤

−− −−

⎟

⎜

⎣

⎦⎣ ⎦

⎟

⎜

⎝⎠

⎟

∑

∑∑

∼

(4.3)

Portanto, neste caso, utiliza-se o amostrador de Gibbs, pois a

distribuição (4.3) possui uma forma fechada, ou seja, a distribuição a posteriori

condicional (4.3) é conhecida.

Agora, a posteriori condicional completa para o parâmetro

do modelo

(2.1) segundo a sua posteriori conjunta (3.5) com base na equação (2.7), é:

()

(

)

()

,, |,

|,,,

,, |,

pN k yH

pk N yH

Nk yHdk

∫

() ()

{}

|,,, exp 1exp

pk N yH y N kt

⎧⎫

⎪⎪

⎡

∝− −−−

⎨⎬

⎣

⎪⎪

⎩⎭

∑

⎤

⎦

(4.4)

A distribuição (4.4) não possui uma forma fechada, então, para o

parâmetro

, não é possível chegar a uma distribuição posteriori condicional

conhecida e, portanto, é necessário se utilizar o algoritmo de Metropolis-

Hastings para se amostrar desta distribuição.

Analogamente, baseando-se na equação (2.7), a condicional completa

para o parâmetro

do modelo (2.1), segundo a sua posteriori conjunta (3.5), é:

()

(

)

()

,, |,

|,,,

,, |,

pN k yH

pNkyH

Nk yHd

ττ

∫

()

{}

|,,,

exp 2 1 exp .

ei i

pNkyH

yN kt

ττ

+ −

∝

⎧⎫

⎡

⎤

⎪⎪

⎡⎤

⎢

⎥

− + −−−

⎨⎬

⎣⎦

⎢

⎥

⎪⎪

⎣

⎦

⎪⎪

⎩⎭

∑

(4.5)

Por inspeção e comparação com a distribuição Gama, tem-se:

αα=+

()

{}

1exp

ee i i

yN ktββ

⎡

⎤

=+ −−−

⎣

⎦

∑

Logo, a distribuição a posteriori condicional completa para

é uma

Gama, ou seja,

() ()

{}

0 0

|,,, , 1exp

eeeeei

NkyH G y N ktτααββ

⎛⎞

⎟

⎜

⎡

⎤

⎟=+ = + −−−

⎜

⎟

⎣

⎦

⎜

⎟

⎜

⎝⎠

∑

∼

. (4.6)

Portanto, neste caso, utiliza-se o amostrador de Gibbs, pois a

distribuição (4.6) possui uma forma fechada, ou seja, a distribuição a posteriori

condicional (4.6) é conhecida.

4.2 Obtenção das posterioris condicionais completas para o modelo de

Marion

Utilizando a mesma linha de raciocínio apresentada anteriormente para o

modelo (2.1), a posteriori condicional completa para o parâmetro

do modelo

(2.2), segundo a sua posteriori conjunta (3.9), baseado na equação (2.7), é dada

por:

()

,,, |,

|,, ,,

,,, |,

pN kb yH

Nkb yHdN

∫

()

{

}

|,, ,, exp 1 exp .

pN kb yH y N kt

⎧⎫

⎪⎪

⎡

⎤

∝− −−−

⎨⎬

⎢

⎥

⎣

⎦

⎪⎪

⎩⎭

∑

(4.7)

Considerando-se

e como vetor dos elementos

para o tempo i, com

1 exp( )= −

1, 2, ,=…in

{

}

,,,= …

yyy y, pode-se definir (4.7)

como:

()()()

[]

() () ()

000

0000

000

|,, ,, exp

exp ' ' ' '

exp ' ' ' ' '

pN kb yH y Nc y Nc

yy yN c cNy cN Nc

yy yN c cNy Ncc

cc cy cc N cc cy

−−

⎧⎫

⎪⎪

∝−− −

⎨⎬

⎪⎪

⎩⎭

⎧⎫

⎪⎪

∝− − − +

⎨⎬

⎪⎪

⎩⎭

⎧⎫

⎪⎪

⎡⎤

∝− − − +

⎨⎬

⎢⎥

⎣⎦

⎪⎪

⎩⎭

⎧⎫

⎪⎪

⎡

⎤⎡

∝− − −

⎨⎬

⎤

⎢

⎥⎢ ⎥

⎣

⎦⎣

⎪

⎩

⎪

⎭

⎦

)

−

que é o núcleo de uma distribuição normal. Portanto, a posteriori condicional

completa para

é uma normal, ou seja,

( ) () ()

(

|,, ,, ' ', '

N kb yH N cc cy ccτ

−

∼ , (4.8)

em que

()

'1exp

⎡⎤

= −−

⎢⎥

⎣⎦

∑

cc kt

()

1exp

kt y

cc cy

−

⎡

⎤

−−

⎢

⎥

⎣

⎦

⎡

⎤

−−

⎢

⎥

⎣

⎦

∑

Assim, a distribuição posteriori condicional completa para

do modelo (2.2)

é uma normal, ou seja,

()

1exp

|,, ,, ,

1exp 1exp

kt y

Nkb yH N

kt kt

⎛⎞

⎡⎤

⎟

⎜

−−

⎟

⎜

⎢⎥

⎟⎣⎦

⎜

⎟

⎜

⎟

⎜

⎟

⎜

⎟

⎡

⎤⎡ ⎤

⎟

⎜

−− −−

⎟

⎜

⎢

⎥⎢ ⎥

⎟⎟

⎣

⎦⎣ ⎦

⎜

⎝⎠

⎟

∑

∑∑

∼

. (4.9)

Portanto, neste caso utiliza-se o amostrador de Gibbs, pois a distribuição

(4.9) possui uma forma fechada, ou seja, a distribuição a posteriori condicional

(4.9) é conhecida.

As posterioris condicionais completas dos parâmetros

do modelo

(2.2), segundo a sua posteriori conjunta (3.9) por meio da equação (2.7), são,

respectivamente:

e kb

()

,,, |,

|,,,,

,,, |,

pN kb yH

pk N b yH

Nkb yHdk

∫

()

{}

| , , , , exp 1 exp ,

pk N b yH y N kt

⎧⎫

⎪⎪

⎡⎤

∝− −−−

⎨⎬

⎢⎥

⎣⎦

⎪⎪

⎩⎭

∑

(4.10)

()

,,, |,

|,,,,

,,, |,

pN kb yH

pbN k yH

Nkb yHdb

∫

()

{

}

|,,,, exp 1exp

pbN k yH y N kt

⎧⎫

⎪⎪

⎡

⎤

∝− −−−

⎨⎬

⎢

⎥

⎣

⎦

⎪⎪

⎩⎭

∑

(4.11)

As distribuições (4.10) e (4.11) não possuem formas fechadas então não

é possível chegar a uma distribuição condicional conhecida para os parâmetros

e, portanto, é necessário utilizar o algoritmo de Metropolis-Hastings para

se amostrar destas distribuições.

e kb

Novamente, por intermédio da equação (2.7), a posteriori condicional

completa para o parâmetro

do modelo (2.2), segundo a sua posteriori

conjunta (3.9), é:

()

(

)

()

,,, |,

|,,,,

,,, |,

pN kb yH

pNkbyH

Nkb yHd

ττ

∫

()

{}

|,,,,

exp 2 1 exp .

ei i

pNkbyH

yN kt

ττ

+ −

∝

⎧⎫

⎡

⎤

⎪⎪

⎡⎤

⎢

⎥

− + −−−

⎨⎬

⎢⎥

⎣⎦

⎢

⎥

⎪⎪

⎣

⎦

⎪⎪

⎩⎭

∑

(4.12)

Por inspeção e comparação com a distribuição Gama, tem-se:

αα=+

()

{

}

1exp

ee i i

yN ktββ

⎡

⎤

=+ −−−

⎢

⎥

⎣

⎦

∑

Sendo assim, a distribuição posteriori condicional completa para

uma Gama, ou seja,

()

{}

** **

0 0

|,,,, , 1exp

eeeeei

NkbyH G y N ktτααββ

⎛⎞

⎟

⎡

⎤

⎜

⎟

=+ = + −−−

⎜

⎟

⎢

⎥

⎜

⎣

⎦

⎟

⎜

⎝⎠

∑

∼

.(4.13)

Portanto, neste caso, utiliza-se o amostrador de Gibbs, pois a

distribuição (4.13) possui uma forma fechada, ou seja, a distribuição a posteriori

condicional (4.13) é conhecida.

4.3 Obtenção das posterioris condicionais completas para o modelo de

Cabrera

Por meio da equação (2.7), a posteriori condicional completa para o

parâmetro

do modelo (2.3), segundo a sua posteriori conjunta (3.13), é dada

por:

()

(

)

()

110

110 1

,,, |,

|,,,,

,,, |,

pN kk yH

Nkk yHdN

∫

() ()

{}

110 1 1 0

|,,,, exp 1exp .

eii

pN kk yH y N kt kt

⎧⎫

⎪⎪

⎡⎤

∝− −−−−

⎨⎬

⎣⎦

⎪⎪

⎩⎭

∑

(4.14)

Considerando-se

{

}

,,,= …

yyy y, ,

= −

yykt

como vetor dos

elementos

, e como vetor dos elementos para o

tempo i, com

pode-se definir (4.14) como:

1 exp( )= −

tdkd

1, 2, ,=…i

()

()()

110 1 1

1111

111

|,,,, exp

'' ''

exp

pN kk yH y Nd y Nd

yy yNd dNy dNNd

yy yNd dNy Ndd

⎧⎫

⎪⎪

∝− − −

⎨⎬

⎪⎪

⎩⎭

⎧⎫

⎪⎪

⎡

⎤

∝− − − +

⎨⎬

⎢

⎥

⎣

⎦

⎪⎪

⎩⎭

⎧⎫

⎪⎪

⎡⎤

∝− − − +

⎨⎬

⎢⎥

⎣⎦

⎪⎪

⎩⎭

() () () ()

110 1 1

''' ''

|,,,, exp

Nkk yH N dd dy ddN dd dy

−−

⎧⎫

⎪⎪

⎡

⎤⎡

∝− − −

⎨⎬

⎤

⎢

⎥⎢ ⎥

⎣

⎦⎣

⎪⎪

⎩⎭

⎦

)

−

que é o núcleo de uma distribuição normal. Portanto, a posteriori condicional

completa para

é uma normal, ou seja,

( ) () ()

(

110

|,,,, ,

Nkk yH Ndd dy ddτ

−

∼ , (4.15)

em que

()

1exp

⎡⎤

= −−

⎣⎦

∑

dd kt

()

(

)

(

)

()

1exp

ii i

kt y kt

−

⎡

⎤

−− −

⎣⎦

⎡⎤

−−

⎣⎦

∑

De acordo com (4.15), a distribuição posteriori condicional completa

para

do modelo (2.3) pode ser reescrita como:

()

()( )

()

110

1exp

|,,,, ,

1exp 1exp

ii i

kt y kt

Nkk yH N

kt kt

⎛⎞

⎡⎤

−− −

⎜

⎟

⎣⎦

⎜

⎟

⎜

⎟

⎜

⎟

⎜

⎟

⎜

⎟

⎡

⎤⎡

−− −−

⎟⎜

⎤

⎣

⎦⎣

⎟

⎜

⎟

⎜

⎝⎠⎟

⎜

∑

∑∑

∼

⎦

. (4.16)

Portanto, neste caso, utiliza-se o amostrador de Gibbs, pois a

distribuição (4.16) possui uma forma fechada, ou seja, a distribuição a posteriori

condicional (4.16) é conhecida.

Baseando-se, novamente, na equação (2.7), a posteriori condicional

completa para o parâmetro

do modelo (2.3), segundo a sua posteriori

conjunta (3.13), é dada por:

()

(

)

()

110

110 1

,,, |,

|,,,,

,,, |,

pNkk yH

pk N k yH

Nkk yHdk

∫

() ()

{}

110 1 1 0

|,,,, exp 1exp .

(4.17)

eiii

pk N k yH y N kt kt

⎧⎫

⎪⎪

⎡⎤

∝− −−−−

⎨⎬

⎣⎦

⎪⎪

⎩⎭

∑

A distribuição (4.17) não possui uma forma fechada, então, para o

parâmetro

não é possível chegar a uma distribuição condicional conhecida e,

portanto, é necessário se utilizar o algoritmo de Metropolis-Hastings para se

amostrar desta distribuição.

Analogamente, com base na equação (2.7), a posteriori condicional

completa para o parâmetro

do modelo (2.3), segundo a sua posteriori

conjunta (3.13), é dada por:

()

110

011

110 0

,,, |,

|,,,,

,,, |,

pN kk yH

pk N k yH

Nkk yHdk

∫

() ()

{}

011 1 1 0

|,,,, exp 1exp .

(4.18)

eiii

pk N k yH y N kt kt

⎧⎫

⎪⎪

⎡⎤

∝− −−−−

⎨⎬

⎣⎦

⎪⎪

⎩⎭

∑

Considerando-se

{

}

,,,= …

yyy y,

()

1exp

⎡

⎤

= −−−

⎣

⎦

Nkt

como vetor dos elementos

e t como um vetor dos tempos, pode-se definir

(4.18) como:

()

()()

()

** **

011 0 0

** ** ** **

00 00

** ** ** **

00 0

|,,,, exp

'' ''

exp

'' '

exp

pk N k yH y kt y kt

yy ykttky tkkt

yy ykttky ktt

ktttyt

−

⎧⎫

⎪⎪

∝− − −

⎨⎬

⎪⎪

⎩⎭

⎧⎫

⎪⎪

⎡

⎤

∝− − − +

⎨⎬

⎢

⎥

⎣

⎦

⎪⎪

⎩⎭

⎧⎫

⎪⎪

⎡⎤

∝− − − +

⎨⎬

⎢⎥

⎣⎦

⎪⎪

⎩⎭

⎡⎤

∝− −

⎢⎥

⎣⎦

() ()

'''

tk tt ty

−

⎧⎫

⎪⎪

⎡

⎤

−

⎨⎬

⎢

⎥

⎣

⎦

⎪⎪

⎩⎭

que é o núcleo de uma distribuição normal. Portanto, a posteriori condicional

completa para

é uma normal, ou seja,

()()

(

011

|, ,,, ,

kkN yH Ntt ty ttτ

−

∼

()

)

−

, (4.19)

em que

∑

tt t

()

(

)

{

}

1exp

ii i

ty N kt

tt ty

−

⎡

⎤

−−−

⎣

⎦

∑

Logo, a distribuição posteriori condicional completa para

pode ser reescrita

como:

()

{

}

011

1exp

|,,,, ,

ii i

ty N kt

kNk yH N

⎛⎞

⎡

⎤

−−−

⎟

⎜

⎣⎦

⎟

⎜

⎟

⎜

⎟

⎜

⎟

⎜

⎟

⎜

⎟

⎜

⎟

⎜

⎟

⎜⎝⎠

∑

∑∑

∼

. (4.20)

Portanto, neste caso utiliza-se o amostrador de Gibbs, pois a distribuição

(4.20) possui uma forma fechada, ou seja, a distribuição a posteriori condicional

(4.20) é conhecida.

Analogamente, com base na equação (2.7), a condicional completa para

o parâmetro

do modelo (2.3), segundo a sua posteriori conjunta (3.13), é:

()

110

,,, |,

|,,,,

,,, |,

pNkk yH

pNkkyH

Nkk yHd

ττ

∫

()

{}

110

|,,,,

exp 2 1 exp .

ei ii

pNkkyH

yN kt kt

ττ

+ −

∝

⎧⎫

⎡⎤

⎪⎪

⎡⎤

⎢⎥

− + −−−−

⎨⎬

⎣⎦

⎢⎥

⎪⎪

⎣⎦

⎪⎪

⎩⎭

∑

(4.21)

Por inspeção e comparação com a distribuição Gama tem-se:

***

αα=+

()

{}

***

1exp

ee i i i

yN kt ktββ

⎡⎤

=+ −−−−

⎣⎦

∑

Logo, a distribuição posteriori condicional completa para

é uma

Gama, ou seja,

() ()

{}

*** ***

110 110

|,,,, , 1exp .

eeeeeii

NkkyH G y N kt ktτααββ

⎛⎞

⎟

⎜

⎡⎤

⎟

=+ = + −−−−

⎜

⎟

⎣⎦

⎜

⎟

⎜

⎝⎠

∑

∼

(4.22)

Portanto, neste caso, utiliza-se o amostrador de Gibbs, pois a

distribuição (4.22) possui uma forma fechada, ou seja, a distribuição a posteriori

condicional (4.22) é conhecida.

4.4 Análise dos dados

No presente trabalho, a monitoração da convergência das cadeias foi

realizada por meio de métodos formais e informais (seção 3.3). A monitoração

formal foi realizada por meio do critério de Geweke e do critério de Raftery &

Lewis. Para o critério de Geweke, existirão evidências contra a convergência da

cadeia, se o p-valor estimado for menor que um nível de significância pré-fixado

pelo pesquisador, enquanto que, no critério de Raftery & Lewis, caracterizam

como não convergência da cadeia se o fator de dependência (FD) apresentar

valor maior que “cinco”. Com base na Tabela 1, pode-se verificar que o critério

de Geweke não apresentou nenhuma evidência contra a convergência dos

parâmetros dos modelos, a 5% e, com relação ao critério de Raftery & Lewis

observa-se que os fatores de dependência para cada parâmetro dos modelos são

inferiores a “cinco”, o que caracteriza a convergência das cadeias. Portanto, os

modelos (2.1). (2.2) e (2.3) convergiram.

TABELA 1 –

Fator de dependência do critério de Raftery & Lewis (FD) e p-

valor do critério de Geweke dos modelos (2.1), (2.2) e (2.3).

Modelo Parâmetro FD Geweke p-valor

1,0109 0,0636

1,0534 0,9899

Stanford & Smith

0,9907 0,9982

1,0109 0,9565

2,5056 0,9925

1,0320 0,9845

Marion

0,9907 0,9834

1,1439 0,8953

0,9704 0,9991

Cabrera

1,0109 0,9552

1,0320 0,9940

A monitoração informal da convergência da cadeia de cada parâmetro

dos modelos (2.1), (2.2) e (2.3) se realizou por meio de técnicas gráficas

baseadas nas trajetórias (Running Mean) e traços das cadeias. Analisando-se as

Figuras 2, 3 e 4, verifica-se que todas as cadeias convergiram, pois, ao observar-

se a trajetória das cadeias ao longo das iterações, verifica-se que os gráficos,

após um período inicial, apresentam, repetidamente, o mesmo comportamento

qualitativo e quantitativo. Nas Figuras 2, 3 e 4, por intermédio do traço, observa-

se que as cadeias se encontram num estado de equilíbrio, ou seja, as cadeias

convergiram. Portanto, os modelos (2.1). (2.2) e (2.3) atingiram a convergência.

Trajetória da cadeia do parâmetro Traço da cadeia do parâmetro

FIGURA 2: Trajetórias e traços das cadeias dos parâmetros

modelo de Stanford & Smith.

Trajetória da cadeia do parâmetro Traço da cadeia do parâmetro

FIGURA 3: Trajetórias e traços das cadeias dos parâmetros

Nkb

modelo de Marion.

Trajetória da cadeia do parâmetro Traço da cadeia do parâmetro

FIGURA 4: Trajetórias e traços das cadeias dos parâmetros

Nk k

11 0

modelo de Cabrera.

Os resultados descritivos das posterioris obtidos no processo de ajuste

dos modelos não lineares estão apresentados na Tabela 2. De acordo com a

Tabela 2, pode-se afirmar que as estimativas dos intervalos HPD mostraram-se

adequadas ao ajuste dos modelos, pois foram obtidas a partir da distribuição

empírica dos parâmetros. Observa-se também que todos os parâmetros dos

modelos (2.1), (2.2) e (2.3) foram significativos, pois seus respectivos intervalos

HPD não continham o valor zero em sua extensão. Pode-se notar, ainda, que os

erros de Monte Carlo apresentaram valores baixos, confirmando a boa qualidade

das estimativas dos parâmetros.

TABELA 2 – Média e desvio padrão a posteriori dos parâmetros (Par.) dos

modelos (2.1), (2.2) e (2.3), com seus respectivos intervalos HPD 95%, erro de

Monte Carlo (Erro MC).

HPD 95%

Modelo Par.

Média

posteriori

LI LS

Desvio

padrão

Erro MC

150,8158 148,2554 153,5245 1,3606 0,01953957

0,1933 0,1806 0,2046 0,0062 0,01581145

Stanford

& Smith

0,0267 0,0093 0,0486 0,0106 0,01581159

149,4726 145,8477 152,8670 1,7807 0,02470511

0,1932 0,1812 0,2052 0,0062 0,01581146

1,0171 1,0000 1,0497 0,0158 0,01581183

Marion

0,0245 0,0075 0,0435 0,0101 0,01581158

256,8945 247,5555 265,7531 4,6515 0,2484605

0,1009 0,0963 0,1057 0,0024 0,0158114

-3,3107 -3,5531 -3,0576 0,1241 0,01584023

Cabrera

0,0446 0,0141 0,0776 0,0173 0,01581184

Os resultados da comparação dos modelos por meio do Fator de Bayes

encontram-se na Tabela 3. Se o Fator de Bayes for maior que 1, o modelo 1 é

mais plausível ou verossímil do que o modelo 2, então, opta-se por escolher o

modelo 1 para predizer o fenômeno em questão. Por intermédio do Fator de

Bayes (Tabela 3), compararam-se os modelos de Stanford & Smith vs Cabrera e

Marion vs Cabrera os resultados foram, 0,0507 e 0,0307, respectivamente. Com

isso, pode-se afirmar que o modelo de Cabrera se mostrou melhor que os demais

no ajuste aos dados. O Fator de Bayes para a comparação Stanford & Smith vs

Marion apresentou resultado de 1,6504, isto é, o modelo de Stanford & Smith é

mais plausível ou verossímil que o modelo de Marion. Portanto, o modelo que

melhor se ajustou aos dados de mineralização de nitrogênio foi o modelo de

Cabrera, seguido pelo modelo de Stanford & Smith e, por último, o de Marion.

TABELA 3 – Fator de Bayes (FB) dos modelos (2.1), (2.2) e (2.3).

Modelos confrontados FB

Stanford & Smith vs Marion 1,6504

Stanford & Smith vs Cabrera 0,0507

Marion vs Cabrera 0,0307

A comparação de modelos, por meio do BIC consiste na escolha do

modelo que apresentar o maior valor de BIC. Portanto, com base nos resultados

apresentados na Tabela 4, o modelo que melhor se ajustou ao conjunto de dados

de mineralização do nitrogênio foi o modelo de Cabrera, seguido pelo modelo de

Stanford & Smith e, por último, o de Marion.

TABELA 4 – Critério de Informação Bayesiano (BIC) dos modelos (2.1), (2.2)

e (2.3).

Modelos BIC

Stanford & Smith -58,77588

Marion -62,17582

Cabrera -55,20933

Nas Figuras 5, 6 e 7 são apresentados, respectivamente, as densidades e

os histogramas dos parâmetros dos modelos (2.1), (2.2) e (2.3). Com base na

simetria das densidades e dos histogramas podem-se propor novas prioris para

estudos futuros. Para os histogramas e densidades simétricas pode-se sugerir a

distribuição Normal ou a distribuição t de Student e, para as assimétricas

poderia-se propor a distribuição Gama.

Densidade do parâmetro Histograma do parâmetro

FIGURA 5: Densidades e histogramas dos parâmetros

do modelo

de Stanford & Smith.

Densidade do parâmetro Histograma do parâmetro

FIGURA 6: Densidades e histogramas dos parâmetros

Nkb

modelo de Marion.

Densidade do parâmetro Histograma do parâmetro

FIGURA 7: Densidades e histogramas dos parâmetros

Nk k

11 0

modelo de Cabrera.

Outro aspecto de grande importância relacionado ao algoritmo de

Metropolis Hastings é a taxa de aceitação, estimada em 0,44 para todos os

modelos que, de acordo com Gelman et al. (2003), deve ser em torno de 0,23 e

0,44, dependendo do número de parâmetros.

Do ponto de vista prático, podem-se predizer as quantidades de N

mineralizadas ao longo do tempo, em solos sob diferentes sistemas de manejo,

com o objetivo de levantar informação sobre o índice de mineralização de solos

representativos de uma região qualquer. Todas essas informações são

importantes para definir com maior precisão a quantidade de adubos

nitrogenados e orgânicos, de modo a obter produtividades maiores das culturas e

aumentar a eficiência de uso de

(nitrato) e (amônia) no solo.

Isso implica na possibilidade de se evitar a contaminação dos mananciais com

resíduos de

e , uma vez que, em excesso (acima de 10 e 0,02

mgL

N-NO

N-NH

N-NO

N-NH

-1

, respectivamente), essas formas de nitrogênio causam danos à saúde dos

homens e dos animais.

A precisão na predição de N potencialmente mineralizável, de fato,

representa uma maior possibilidade de se definir com precisão a dose de N a ser

adicionado ao solo para as culturas. Isso implica em economia, uso racional de

recursos escassos e possibilidade de melhoria na qualidade ambiental.

5 CONCLUSÕES

Diante dos resultados obtidos, pode-se concluir que:

a) Nas condições em que foram realizados os estudos (tipo de solo, região,

presença de calagem, etc.), o modelo que proporcionou melhor qualidade

de ajuste ao conjunto de dados de mineralização de nitrogênio foi o

modelo de Cabrera, tanto utilizando o Fator de Bayes quanto o Critério de

Informação Bayesiano.

b) Nas condições em que foram realizados os estudos (tipo de solo, região,

presença de calagem, etc), o modelo que proporcionou pior qualidade de

ajuste ao conjunto de dados de mineralização de nitrogênio foi o modelo

de Marion, tanto utilizando o Fator de Bayes quanto o Critério de

Informação Bayesiano.

c) A metodologia bayesiana mostrou-se de fácil aplicação e compreensão. As

estatísticas descritivas das distribuições marginais são facilmente obtidas

para cada parâmetro, como, por exemplo, no cálculo de intervalos de

credibilidades.

d) A inferência bayesiana é uma metodologia de grande potencial para o

estudo da mineralização do nitrogênio, pois considera a incerteza existente

sobre todos os parâmetros dos modelos estudados, além de possibilitar a

inclusão de informações por meio do uso de distribuições a priori

informativas.

Em estudos futuros, pretende-se levantar informação sobre o índice de

mineralização do nitrogênio em solos representativos da região do Sul de Minas

Gerais, procurando quantificar os estoques de nitrogênio, a fim de determinar se

há necessidade ou não de adubação nitrogenada no solo. Pode-se assim, avaliar a

qualidade do solo com base no estoque de nitrogênio e na composição da

matéria orgânica do solo. Outro fato a ser verificado nestes estudos seria definir

práticas de manejo ou, melhor, definir doses e ou aporte de nitrogênio em solo

que melhora a qualidade e aumenta o estoque de matéria orgânica em solo e,

consequentemente, avaliar o impacto ambiental do nitrogênio no que diz respeito

à qualidade do solo e suas inter-relações com a matéria orgânica do solo.

No tocante à inferência bayesiana, pretende-se propor novas

distribuições a priori (informativa) para os parâmetros dos modelos presentes na

literatura.

6 REFERÊNCIAS BIBLIOGRÁFICAS

BERKHOF, J.; MECHELEN, I. V.; GELMAN, A. A bayesian approach to the

selection and testing of mixture models.

Statistica Sinica, Taipei, v. 13, n. 2, p.

423-442, Apr. 2003.

BLACK, C. A.

Soil-plant relationships. New York: John Wiley, 1968. 792 p.

BOX, G. E. P.; TIAO, G. C.

Bayesian inference in statistical analysis. New

York: John Wiley & Sons, 1992. 588 p.

BROADBENT, F. E. Empirical modeling of soil nitrogen mineralization.

Soil

Science

, Baltimore, v. 141, n. 3, p. 208-213, Mar. 1986.

CABRERA, M. L. Modeling the flush of nitrogen mineralization caused by

drying and rewetting soils.

Soil Science Society of America Journal, Madison,

v. 57, n. 1, p. 63-66, Jan. /Feb. 1993.

CAMARGO, F. A. de O.; GIANELLO, C.; TEDESCO, M. J.; RIBOLDI, J.;

MEURER, E. J.; BISSANI, C. A. Empirical models to predict soil nitrogen

mineralization.

Ciência Rural, Santa Maria, v. 32, n. 3, p. 393-399, maio/jun.

2002.

CASELLA, G.; GEORGE, E. Explaning the gibbs sampler.

The American

Statitician

, Alexandria, v. 46, n. 3, p. 167-174, Aug. 1992.

CARLIN, B. P.; LOUIS, T. A.

Bayes and empirical bayes methods for data

analysis

. 2. ed. London: Chapman and Hall, 2000. 419 p.

CHEN, M. H.; SHAO, Q. M. Monte Carlo estimation of Bayesian credible and

HPD intervals.

Journal of Computational and Graphical Stattistics,

Alexandria, v. 8, n. 1, p. 69-92, Mar. 1999.

CHEW, W. Y.; WILLIAMS, C. N.; JOSEPH, K. T.; RAMLI, K. Studies on the

availability to plants of soil nitrogen in Malaysian tropical oligotrophic peat. I –

Effect of liming and pH.

Tropical Agriculture, London, v. 53, p. 9-78, Jan.

1976.

CHIB, S.; GREENBERG, E. Understanding the Metropolis-Hastings algorithm.

The American Statistician, Alexandria, v. 49, n. 4, p. 327-335, Nov. 1995.

COWLES, M. K.; CARLIN, B. P. Markov chain Monte Carlo convergence

diagnostics: a comparative review.

Journal of the American Statistical

Association

, Alexandria, v. 91, n. 434, p. 883-904, June 1996.

DRAPER, N. R.; SMITH, H.

Applied regression analysis. 3. ed. New York:

John Wiley e Sons, 1998. 706 p.

EHLERS, R. S.

Introdução à inferência bayesiana. 2. ed. Curitiba, 2003. 68 p.

FERGUSON, R. B.; SHAPIRO, C. A.; HERGERT, G. W.; KRANZ, W. L.;

KLOCKE, N. L.; KRULL, D. H. Nitrogen and irrigation management practices

to minimize nitrate leaching from irrigate corn.

Journal of Production

Agriculture,

Madison, v. 4, n. 2, p. 186-192, Apr./June 1991.

GAMERMAN, D.

Markov Chain Monte Carlo: Stochastic Simulation For

Bayesian Inference. Londres: Chapman & Hall, 1997. 245 p.

GELFAND, A. E.; SMITH, A. F. M. Sampling based approaches for calculating

marginal densities.

Journal of the American Statistical Association, London,

v. 85, n. 410, p. 398-409, June 1990.

GELMAN, A.; CARLIN, J. B.; STERN, H. S.; RUBIN, D. B.

Bayesian data

analysis

. London: Chapman and Hall, 2003. 668 p.

GELMAN, A.; RUBIN, D. B. Inference from iterative simulation using multiple

sequences.

Statistical Science, Hayward, v. 7, n. 4, p. 457-511, Nov. 1992.

GEMAN, S.; GEMAN, D. Stochastic relaxation, Gibbs distributions and the

Bayesian restoration of images.

IEEE Transactions on pattern Analysis and

Machine Intellingence

, Los Alamitos, v. 6, n. 6, p. 721-741, 1984.

GEWEKE, J. Evaluating the accuracy of sampling-based approaches to the

calculation of posterior moments. N: BERNARDO, J. M.; BERGER, J. O.;

DAWID, A. P.; SMITH, A. F. M. (Ed. ).

Bayesian Statistics 4. New York:

Oxford University Press, 1992. 526 p.

GUPTA, U. C.; REUSZER, H. W. Effect of plant species on the aminoacid

content and nitrification of soil organic matter.

Soil Science, Baltimore, v. 104,

n. 6, p. 395-400, 1967.

Han, C.; Carlin, B. MCMC methods for computing Bayes Factors: a

comparative review.

Journal of the American Statistical Association,

Alexandria, v. 96, n. 455, p. 1122-1132, Sept. 2001.

HASTINGS, W. K. Monte carlo sampling methods using Markov chains and

their applications.

Biometrika, London, v. 57, n. 1, p. 97-109, Apr. 1970.

HEIDELBERGER, P.; WELCH, P. Simulation run length control in the

presence of an initial transient.

Operations Research, Landing, v. 31, n. 6, p.

1109-1144, 1983.

INOBUSHI, K.; WADA, H.; TAKAI, Y. Easily decomposable organic matter in

paddy soil. VI. Kinetics of nitrogen mineralization in submerged soils.

Soil

Science and Plant Nutrition,

Tokyo, v. 31, n. 4, p. 563-572, Dec. 1985.

JONES, A. Estimation of an active fraction soil nitrogen.

Communications in

Soil Science and Plant Analysis,

New York, v. 15, n. 1, p. 23-32, 1984.

JUMA, N. G.; PAUL, E. A.; MARY, B. Kinetic analysis of net mineralization in

soil.

Soil Science Society of America Journal, Madison, v. 48, n. 4, p. 465-472,

July/Aug. 1984.

MARION, G. M.; KUMMEROW, J.; MILLER, P. C. Predicting nitrogen

mineralization in chaparral soils.

Soil Science Society of America Journal,

Madison, v. 45, n. 5, p. 956-961, Sept./Oct. 1981.

METROPOLIS, N.; ROSENBLUTH, A. W.; ROSENBLUTH, M. N.; TELLER,

A. H.; TELLER, E. Equations of state calculations by fast computing machines.

Journal of Chemical Phisics, Woodbury, v. 21, n. 6, 1087-1092, 1953.

MOLINA, J. A. E.; CLAPP, C. E.; LARSON, W. E. Potentially mineralizable

nitrogen in soil: the simple exponential model does not apply for the first 12

weeks of incubation.

Soil Science Society of American Journal, Madison, v.

44, n. 2, p. 442-443, Mar./Apr. 1980.

PAULINO, C. D.; TURKMAN, A. A.; MURTEIRA, B.

Estatística Bayesiana,

Lisboa: Fundação Calouste Gulbenkian, 2003. 446 p.

PEREIRA, J. M.

Predição do nitrogênio mineralizado em latossolo pelo uso

de modelos não lineares

. 2004. 82 p. Dissertação (Mestrado em Estatística e

Experimentação Agropecuária) – Universidade Federal de Lavras, Lavras, MG.

PEREIRA, J. M.; MUNIZ, J. A.; SILVA, C. A. Nonlinear models to predict

nitrogen mineralization in a oxisol.

Scientia Agricola, Piracicaba, v. 62, n. 4, p.

395-400, July/Aug. 2005.

RAFTERY, A. L.; LEWIS, S. comment: one long run with diagnostics:

implementation strategies for Markov chain Monte Carlo.

Statistical Science,

Hayward, v. 7, n. 4, 493-497, Nov. 1992.

RAIJ, B. van.

Fertilidade do solo e adubação. Piracicaba: Ceres, 1991. 343 p.

R DEVELOPMENT CORE TEAM.

The R Manuals: Current Version: 1. 9. 1.

2004. Disponível em: <http://www.r-project.org>. Acesso em: 2006.

ROSA, G. M.

Análise Bayesiana de modelos lineares mistos robustos via

Amostrador de Gibbs

. 1998. 57 p. Tese (Doutorado em Estatística e

Experimentação Agronômica) – Universidade de São Paulo, Piracicaba, SP.

SAINJU, U. M.; SINGH, B. P.; RAHMAN, S.; REDDY, V. R. Soil nitrat-

nitrogen under tomato following tillage, cover cropping, and nitrogen

fertilization.

Journal of Environmental Quality, Madison, v. 28, n. 6, p. 1837-

1844, Nov./dez. 1999.

SILVA, C. A.; VALE, F. R.; GUILHERME, L. R. Efeito da calagem na

mineralização do nitrogênio em solos de Minas Gerais.

Revista Brasileira

Ciência do Solo

, Campinas, v. 18, n. 3, p. 471-476, set./dez. 1994.

SORENSEN, D.

Gibbs sampling in quantitative genetics. Foulun, 1996. 186

STANFORD, G.; SMITH, S. J. Nitrogen mineralization potentials of soil.

Soil

Science of America Journal

, Madison, v. 36, n. 3, p. 465-472, May/June 1972.

STANFORD, G.; SMITH, S. J. Estimating potentially mineralizable soil

nitrogen from a chemical index of soil nitrogen availability.

Soil Science,

Baltimore, v. 122, n. 2, p. 71-76, 1976.

UNLU, K.; OZENIRLER, G.; YURTERI, C. Nitrogen fertilizer leaching from

cropped and irrigated sandy soil in Central Turkey.

European Journal of Soil

Science

, Oxford, v. 50, n. 4, p. 609-620, Dec. 1999.

ZEGER, S. L.; KARIM, M. R. Generalized linear models with random effects; a

Gibbs sampling approach.

Journal of the American Statistical Association,

Alexandria, v. 86, n. 413, p. 79-86, Mar. 1991.

7 ANEXO

PROGRAMAS Página

PROGRAMA 1

Rotina R para gerar as cadeias via MCMC

do modelo de Stanford & Smith................... 57

PROGRAMA 2

Rotina R para realizar o “burn-in” e o

“thin” do modelo de Stanford & Smith........ 61

PROGRAMA 3

Rotina R para gerar as cadeias via MCMC

do modelo de Marion.................................... 63

PROGRAMA 4

Rotina R para realizar o “burn-in” e o

“thin” do modelo de Marion......................... 69

PROGRAMA 5

Rotina R para gerar as cadeias via MCMC

do modelo de Cabrera...................................

PROGRAMA 6

Rotina R para realizar o “burn-in” e o

“thin” do modelo de Cabrera........................ 77

Programa 1: Rotina R para gerar as cadeias via MCMC do modelo de Stanford

& Smith.

dados<-read.table("arquivo de dados", h=T)

y<-dados$y

t<-dados$t

n<-length(y)

# Definindo valores para os Hiperparâmetros #

alpha=1.5

beta=1

# Definindo contador para a taxa de aceitação #

contk=0

taxak<-0

# Assumindo Prioris #

# teta <- cte #

# tau <- Gama(Alpha e Beta) #

# Definindo o número de iterações #

nit<-450000

No<-matrix(0,nit,1)

k<-matrix(0,nit,1)

tau<-matrix(0,nit,1)

# Pontos iniciais (Chutes) #

No[1,1]<-158

k[1,1]<-0.2

tau[1,1]<-4

ir<-0

for (ir in 2:nit)

{

# Condicionais completas #

# Condicional completa para o parâmetro

somaa1<-0

somaa<-0

somab<-0

somaa2<-0

i<-0

for (i in 1:n)

{

somaa=y[i]*(1-exp(-k[ir-1,1]*t[i]))

somaa1=somaa1+somaa

somab=(1-exp(-k[ir-1,1]*t[i]))^2

somaa2=somaa2+somab

}

mediaNo<-somaa1/somaa2

desvNo<-sqrt(tau[ir-1,1]/somaa2)

No1<-rnorm(1, mean=mediaNo, sd=desvNo)

No[ir,1]<-No1

# Condicional completa para o parâmetro

somac<-0

somaa3<-0

i<-0

for (i in 1:n)

{

somac=(y[i]-No[ir,1]*(1-exp(-k[ir-1,1]*t[i])))^2

somaa3=somaa3+somac

}

alpha1=(n/2)+alpha

beta1=(beta + 0.5*somaa3)

tau1<-rgamma(1, shape = alpha1, scale = 1/beta1)

tau[ir,1]<-tau1

# Metropolis para o parâmetro k #

iii<-0

while (iii<3)

{

knew<-rnorm(1,0.2,0.006)

if(abs(knew)<1) iii<-3 else iii<-1

}

somacnew=0

soma3new=0

ii<-0

for (ii in 1:n)

{

somacnew=(y[ii]-No[ir,1]*(1-exp(-knew*t[ii])))^2

soma3new=soma3new+somacnew

}

somacold=0

soma3old=0

i<-0

for (i in 1:n)

{

somacold=(y[i]-No[ir,1]*(1-exp(-k[ir-1,1]*t[i])))^2

soma3old=soma3old+somacold

}

knew1<-exp(-(tau[ir,1]/2)*soma3new)

kold1<-exp(-(tau[ir,1]/2)*soma3old)

if(k[ir-1,1]!=0) p<-min(1,(knew1/kold1)) else p<-1

u<-runif(1)

if(u<=p) k1<-knew else k1<-k[ir-1,1]

k[ir,1]<-k1

if(u<=p) contk<-contk+1

}

# Fim de Amostragem #

# Determinando a taxa de aceitação de k #

taxak<-contk/nit

taxak

# Salvando cadeia de cada Parâmetro #

write(k, "saida_k_Stanford_Smith_450000.txt", ncolumns =1)

write(No, "saida_No_Stanford_Smith_450000.txt", ncolumns =1)

write(tau, "saida_tau_Stanford_Smith_450000.txt", ncolumns =1)

Programa 2: Rotina R para realizar o “burn-in” e o “thin” do modelo de

Stanford & Smith.

No<-read.table("arquivo de dados", h=T)

k<-read.table("arquivo de dados", h=T)

tau<-read.table("arquivo de dados", h=T)

No<-as.matrix(No)

k<-as.matrix(k)

tau<-as.matrix(tau)

# Definindo o “burn-in” e o “thin” #

nit<-length(No)

burn<-50000 # descartando as 50.000 primeiras iterações #

pulo<-100 # salto ou pulo a cada 100 iterações #

# Defina os parâmetros #

No<-No[(burn+1):nit,1]

k<-k[(burn+1):nit,1]

tau<-tau[(burn+1):nit,1] # Fim do “burn in” #

# Início do Pulo #

nit_novo<-((nit-burn)%/%pulo)

No_pulo<-matrix(0,nit_novo,1)

k_pulo<-matrix(0,nit_novo,1)

tau_pulo<-matrix(0,nit_novo,1)

# Realizando o “thin” (pulo) individualmente para cada parâmetro #

i<-0

irr<-0

for (irr in 1:(nit-burn))

{

if(irr-(irr%/%pulo)*pulo==0) (i<-i+1)

if(irr-(irr%/%pulo)*pulo==0) (No_pulo[i,1]<-No[irr])

}

i<-0

irr<-0

for (irr in 1:(nit-burn))

{

if(irr-(irr%/%pulo)*pulo==0) (i<-i+1)

if(irr-(irr%/%pulo)*pulo==0) (k_pulo[i,1]<-k[irr])

}

i<-0

irr<-0

for (irr in 1:(nit-burn))

{

if(irr-(irr%/%pulo)*pulo==0) (i<-i+1)

if(irr-(irr%/%pulo)*pulo==0) (tau_pulo[i,1]<-tau[irr])

} # Fim do pulo #

# Salvando cadeia de cadeia parâmetro #

write(No_pulo,“saida_No_Stanford_Smith_450000_50000_100_final_4000.txt"

, ncolumns =1)

write(k_pulo, "saida_k_Stanford_Smith_450000_50000_100_final_4000.txt",

ncolumns =1)

write(tau_pulo,“saida_tau_Stanford_Smith_450000_50000_100_final_4000.txt"

, ncolumns =1)

Programa 3: Rotina R para gerar as cadeias via MCMC do modelo de Marion.

dados<-read.table("arquivo de dados", h=T)

y<-dados$y

t<-dados$t

n<-length(y)

# Definindo valores para os Hiperparâmetros #

alpha=1.5

beta=1

# Definindo contador para as taxas de aceitação dos parâmetros k e b #

contk=0

contb=0

taxak<-0

taxab<-0

# Assumindo Prioris #

# teta <- cte #

# tau <- Gama(alpha e Beta) #

# Definindo o número de iterações #

nit<-450000

No<-matrix(0,nit,1)

k<-matrix(0,nit,1)

b<-matrix(0,nit,1)

tau<-matrix(0,nit,1)

# Pontos iniciais (Chutes) #

No[1,1]<-140

k[1,1]<-0.2

b[1,1]<-1

tau[1,1]<-4

ir<-0

for (ir in 2:nit)

{

# Condicionais completas #

# Condicional completa para o parâmetro

somaa1<-0

somaa<-0

somab<-0

somaa2<-0

iv<-0

for (iv in 1:n)

{

somaa=y[iv]*(1-exp(-k[ir-1,1]*(t[iv]^b[ir-1,1])))

somaa1=somaa1+somaa

somab=(1-exp(-k[ir-1,1]*(t[iv]^b[ir-1,1])))^2

somaa2=somaa2+somab

}

mediaNo<-somaa1/somaa2

desvNo<-sqrt(tau[ir-1,1]/somaa2)

No1<-rnorm(1, mean=mediaNo, sd=desvNo)

No[ir,1]<-No1

# Condicional completa para o parâmetro

somac<-0

somaa3<-0

v<-0

for (v in 1:n)

{

somac=(y[v]-No[ir,1]*(1-exp(-k[ir-1,1]*(t[v]^b[ir-1,1]))))^2

somaa3=somaa3+somac

}

alpha1=(n/2)+alpha

beta1=(beta + 0.5*somaa3)

tau1<-rgamma(1, shape = alpha1, scale = 1/beta1)

tau[ir,1]<-tau1

# Metropolis para o parâmetro k #

iii<-0

while (iii<3)

{

knew<-rnorm(1,0.2,0.006)

if(abs(knew)<1) iii<-3 else iii<-1

}

somacnew=0

soma3new=0

ii<-0

for (ii in 1:n)

{

somacnew=(y[ii]-No[ir,1]*(1-exp(-knew*(t[ii]^b[ir-1,1]))))^2

soma3new=soma3new+somacnew

}

somacold=0

soma3old=0

i<-0

for (i in 1:n)

{

somacold=(y[i]-No[ir,1]*(1-exp(-k[ir-1,1]*(t[i]^b[ir-1,1]))))^2

soma3old=soma3old+somacold

}

knew1<-exp(-(tau[ir,1]/2)*soma3new)

kold1<-exp(-(tau[ir,1]/2)*soma3old)

if(k[ir-1,1]!=0) p<-min(1,(knew1/kold1)) else p<-1

u<-runif(1)

if(u<=p) k1<-knew else k1<-k[ir-1,1]

k[ir,1]<-k1

if(u<=p) contk<-contk+1

# Metropolis para o parâmetro b #

jjj<-0

while (jjj<3)

{

bnew<-rnorm(1,1,0.06)

if(bnew>1) jjj<-3 else jjj<-1

}

somadnew=0

soma4new=0

jj<-0

for (jj in 1:n)

{

somadnew=(y[jj]-No[ir,1]*(1-exp(-k[ir,1]*(t[jj]^bnew))))^2

soma4new=soma4new+somadnew

}

somadold=0

soma4old=0

j<-0

for (j in 1:n)

{

somadold=(y[j]-No[ir,1]*(1-exp(-k[ir,1]*(t[j]^b[ir-1,1]))))^2

soma4old=soma4old+somadold

}

bnew1<-exp(-(tau[ir,1]/2)*soma4new)

bold1<-exp(-(tau[ir,1]/2)*soma4old)

if(b[ir-1,1]!=0) p<-min(1,(bnew1/bold1)) else p<-1

u<-runif(1)

if(u<=p) b1<-bnew else b1<-b[ir-1,1]

b[ir,1]<-b1

if(u<=p) contb<-contb+1

}

# Fim de Amostragem #

# Determinando as taxas de aceitação dos parâmetros k e b #

taxak<-contk/nit

taxab<-contb/nit

taxak

taxab

# Salvando cadeia de cada Parâmetro #

write(b, "saida_b_Marion_450000.txt", ncolumns =1)

write(k, "saida_k_Marion_450000.txt", ncolumns =1)

write(No, "saida_No_Marion_450000.txt", ncolumns =1)

write(tau, "saida_tau_Marion_450000.txt", ncolumns =1)

Programa 4: Rotina R para realizar o “burn-in” e o “thin” do modelo de

Marion.

No<-read.table("arquivo de dados", h=T)

k<-read.table("arquivo de dados", h=T)

b<-read.table("arquivo de dados", h=T)

tau<-read.table("arquivo de dados", h=T)

No<-as.matrix(No)

k<-as.matrix(k)

b<-as.matrix(b)

tau<-as.matrix(tau)

# Definindo o “burn-in” e o “thin” #

nit<-length(No)

burn<-50000 # descartando as 50.000 primeiras iterações #

pulo<-100 # salto ou pulo a cada 100 iterações #

# Defina os parâmetros #

No<-No[(burn+1):nit,1]

k<-k[(burn+1):nit,1]

b<-b[(burn+1):nit,1]

tau<-tau[(burn+1):nit,1] # Fim do “burn-in” #

# Início do Pulo #

nit_novo<-((nit-burn)%/%pulo)

No_pulo<-matrix(0,nit_novo,1)

k_pulo<-matrix(0,nit_novo,1)

b_pulo<-matrix(0,nit_novo,1)

tau_pulo<-matrix(0,nit_novo,1)

# Realizando o “thin” (pulo) individualmente para cada parâmetro #

i<-0

irr<-0

for (irr in 1:(nit-burn))

{

if(irr-(irr%/%pulo)*pulo==0) (i<-i+1)

if(irr-(irr%/%pulo)*pulo==0) (No_pulo[i,1]<-No[irr])

}

i<-0

irr<-0

for (irr in 1:(nit-burn))

{

if(irr-(irr%/%pulo)*pulo==0) (i<-i+1)

if(irr-(irr%/%pulo)*pulo==0) (k_pulo[i,1]<-k[irr])

}

i<-0

irr<-0

for (irr in 1:(nit-burn))

{

if(irr-(irr%/%pulo)*pulo==0) (i<-i+1)

if(irr-(irr%/%pulo)*pulo==0) (b_pulo[i,1]<-b[irr])

}

i<-0

irr<-0

for (irr in 1:(nit-burn))

{

if(irr-(irr%/%pulo)*pulo==0) (i<-i+1)

if(irr-(irr%/%pulo)*pulo==0) (tau_pulo[i,1]<-tau[irr])

} # Fim do pulo #

# Salvando cadeia de cadeia parâmetro #

write(No_pulo, "saida_No_Marion_450000_50000_100_final_4000.txt",

ncolumns =1)

write(k_pulo, "saida_k_Marion_450000_50000_100_final_4000.txt", ncolumns

=1)

write(b_pulo, "saida_b_Marion_450000_50000_100_final_4000.txt", ncolumns

=1)

write(tau_pulo, "saida_tau_Marion_450000_50000_100_final_4000.txt",

ncolumns =1)

Programa 5: Rotina R para gerar as cadeias via MCMC do modelo de Cabrera.

dados<-read.table("arquivo de dados", h=T)

y<-dados$y

t<-dados$t

n<-length(y)

# Definindo valores para os Hiperparâmetros #

alpha=1.5

beta=1

# Definindo contador para as taxas de aceitação dos parâmetros k e b #

contk1=0

taxak1<-0

# Assumindo Prioris #

# teta <- cte #

# tau <- Gama(alpha e beta) #

# Definindo o número de iterações #

nit<-450000

N1<-matrix(0,nit,1)

k1<-matrix(0,nit,1)

ko<-matrix(0,nit,1)

tau<-matrix(0,nit,1)

# Pontos iniciais (Chutes) #

N1[1,1]<-260

k1[1,1]<-0.1

ko[1,1]<--2

tau[1,1]<-4

ir<-0

for (ir in 2:nit)

{

# Condicionais completas #

# Condicional completa para o parâmetro

somaa1<-0

somaa<-0

somab<-0

somaa2<-0

iv<-0

for (iv in 1:n)

{

somaa=(y[iv]-ko[ir-1,1]*t[iv])*(1-exp(-k1[ir-1,1]*t[iv]))

somaa1=somaa1+somaa

somab=(1-exp(-k1[ir-1,1]*t[iv]))^2

somaa2=somaa2+somab

}

mediaN1<-somaa1/somaa2

desvN1<-sqrt(tau[ir-1,1]/somaa2)

N11<-rnorm(1, mean=mediaN1, sd=desvN1)

N1[ir,1]<-N11

# Condicional completa para o parâmetro

somac<-0

somaa3<-0

v<-0

for (v in 1:n)

{

somac=(y[v]-N1[ir,1]*(1-exp(-k1[ir-1,1]*t[v]))-ko[ir-1,1]*t[v])^2

somaa3=somaa3+somac

}

alpha1=(n/2)+alpha

beta1=(beta + 0.5*somaa3)

tau1<-rgamma(1, shape = alpha1, scale = 1/beta1)

tau[ir,1]<-tau1

#Condicional completa para o parâmetro

somad<-0

somaa4<-0

somae<-0

somaa5<-0

v<-0

for (v in 1:n)

{

somad=t[v]*(y[v]-N1[ir,1]*(1-exp(-k1[ir-1,1]*t[v])))

somaa4=somaa4+somad

somae=t[v]^2

somaa5=somaa5+somae

}

mediako<-somaa4/somaa5

desvko<-sqrt(tau[ir,1]/somaa5)

ko1<-rnorm(1, mean=mediako, sd=desvko)

ko[ir,1]<-ko1

# Metropolis para o parâmetro

iii<-0

while (iii<3)

{

k1new<-rnorm(1,0.1,0.003)

if(abs(k1new)<1) iii<-3 else iii<-1

}

somacnew=0

soma3new=0

ii<-0

for (ii in 1:n)

{

somacnew=(y[ii]-N1[ir,1]*(1-exp(-k1new*t[ii]))-ko[ir,1]*t[ii])^2

soma3new=soma3new+somacnew

}

somacold=0

soma3old=0

i<-0

for (i in 1:n)

{

somacold=(y[i]-N1[ir,1]*(1-exp(-k1[ir-1,1]*t[i]))-ko[ir,1]*t[i])^2

soma3old=soma3old+somacold

}

k1new1<-exp(-(tau[ir,1]/2)*soma3new)

k1old1<-exp(-(tau[ir,1]/2)*soma3old)

if(k1[ir-1,1]!=0) p<-min(1,(k1new1/k1old1)) else p<-1

u<-runif(1)

if(u<=p) k1[ir,1]<-k1new else k1[ir,1]<-k1[ir-1,1]

if(u<=p) contk1<-contk1+1

}

# Fim de Amostragem #

# Determinando as taxas de aceitação dos parâmetros k e b #

taxak1<-contk1/nit

taxak1

# Salvando cadeia de cadeia parâmetro #

write(N1, "saida_N1_Cabrera_450000.txt", ncolumns =1)

write(k1, "saida_k1_Cabrera_450000.txt", ncolumns =1)

write(ko, "saida_ko_Cabrera_450000.txt", ncolumns =1)

write(tau, "saida_tau_Cabrera_450000.txt", ncolumns =1)

Programa 6: Rotina R para realizar o “burn-in” e o “thin” do modelo de

Cabrera.

N1<-read.table("arquivo de dados", h=T)

k1<-read.table("arquivo de dados", h=T)

ko<-read.table("arquivo de dados", h=T)

tau<-read.table("arquivo de dados", h=T)

N1<-as.matrix(N1)

k1<-as.matrix(k1)

ko<-as.matrix(ko)

tau<-as.matrix(tau)

# Definindo o “burn-in” e o “thin” #

nit<-length(N1)

burn<-50000 # descartando as 50.000 primeiras iterações #

pulo<-100 # salto ou pulo a cada 100 iterações #

# Defina os parâmetros #

N1<-N1[(burn+1):nit,1]

k1<-k1[(burn+1):nit,1]

ko<-ko[(burn+1):nit,1]

tau<-tau[(burn+1):nit,1] # Fim do “burn-in” #

# Início do Pulo #

nit_novo<-((nit-burn)%/%pulo)

N1_pulo<-matrix(0,nit_novo,1)

k1_pulo<-matrix(0,nit_novo,1)

ko_pulo<-matrix(0,nit_novo,1)

tau_pulo<-matrix(0,nit_novo,1)

# Realizando o “thin” (pulo) individualmente para cada parâmetro #

i<-0

irr<-0

for (irr in 1:(nit-burn))

{

if(irr-(irr%/%pulo)*pulo==0) (i<-i+1)

if(irr-(irr%/%pulo)*pulo==0) (N1_pulo[i,1]<-N1[irr])

}

i<-0

irr<-0

for (irr in 1:(nit-burn))

{

if(irr-(irr%/%pulo)*pulo==0) (i<-i+1)

if(irr-(irr%/%pulo)*pulo==0) (k1_pulo[i,1]<-k1[irr])

}

i<-0

irr<-0

for (irr in 1:(nit-burn))

{

if(irr-(irr%/%pulo)*pulo==0) (i<-i+1)

if(irr-(irr%/%pulo)*pulo==0) (ko_pulo[i,1]<-ko[irr])

}

i<-0

irr<-0

for (irr in 1:(nit-burn))

{

if(irr-(irr%/%pulo)*pulo==0) (i<-i+1)

if(irr-(irr%/%pulo)*pulo==0) (tau_pulo[i,1]<-tau[irr])

} # Fim do pulo #

# Salvando cadeia de cadeia parâmetro #

write(N1_pulo, "saida_N1_Cabrera_450000_50000_100_final_4000.txt",

ncolumns =1)

write(k1_pulo, "saida_k1_Cabrera_450000_50000_100_final_4000.txt",

ncolumns =1)

write(ko_pulo, "saida_ko_Cabrera_450000_50000_100_final_4000.txt",

ncolumns =1)

write(tau_pulo, "saida_tau_Cabrera_450000_50000_100_final_4000.txt",

ncolumns =1)

Livros Grátis
( http://www.livrosgratis.com.br )
 
Milhares de Livros para Download:
 
Baixar livros de Administração
Baixar livros de Agronomia
Baixar livros de Arquitetura
Baixar livros de Artes
Baixar livros de Astronomia
Baixar livros de Biologia Geral
Baixar livros de Ciência da Computação
Baixar livros de Ciência da Informação
Baixar livros de Ciência Política
Baixar livros de Ciências da Saúde
Baixar livros de Comunicação
Baixar livros do Conselho Nacional de Educação - CNE
Baixar livros de Defesa civil
Baixar livros de Direito
Baixar livros de Direitos humanos
Baixar livros de Economia
Baixar livros de Economia Doméstica
Baixar livros de Educação
Baixar livros de Educação - Trânsito
Baixar livros de Educação Física
Baixar livros de Engenharia Aeroespacial
Baixar livros de Farmácia
Baixar livros de Filosofia
Baixar livros de Física
Baixar livros de Geociências
Baixar livros de Geografia
Baixar livros de História
Baixar livros de Línguas

Baixar livros de Literatura
Baixar livros de Literatura de Cordel
Baixar livros de Literatura Infantil
Baixar livros de Matemática
Baixar livros de Medicina
Baixar livros de Medicina Veterinária
Baixar livros de Meio Ambiente
Baixar livros de Meteorologia
Baixar Monografias e TCC
Baixar livros Multidisciplinar
Baixar livros de Música
Baixar livros de Psicologia
Baixar livros de Química
Baixar livros de Saúde Coletiva
Baixar livros de Serviço Social
Baixar livros de Sociologia
Baixar livros de Teologia
Baixar livros de Trabalho
Baixar livros de Turismo
 
 