Download PDF
ads:
ALGORITMOS PARA MODELOS DE LIMIAR
USANDO AS DISTRIBUIÇÕES ACUMULADAS
NORMAL E t DE STUDENT
JOSÉ WALDEMAR DA SILVA
2008
ads:
Livros Grátis
http://www.livrosgratis.com.br
Milhares de livros grátis para download.
JOSÉ WALDEMAR DA SILVA
ALGORITMOS PARA MODELOS DE LIMIAR USANDO
AS DISTRIBUIÇÕES ACUMULADAS NORMAL E t DE
STUDENT
Tese apresentada à Universidade Federal de Lavras
como parte das exigências do Programa de Pós-
graduação em Estatística e Experimentação Agro-
pecuária, para a obtenção do título de “Doutor”.
Orientador
Prof. Dr. Júlio Sílvio de Sousa Bueno Filho
LAVRAS
MINAS GERAIS-BRASIL
2008
ads:
Ficha Catalográfica Preparada pela Divisão de Processos Técnicos da
Biblioteca Central da UFLA
Silva, José Waldemar da.
Algoritmos para modelos de limiar usando as distribuições acumuladas
normal e t de Student / José Waldemar da Silva. -- Lavras : UFLA, 2008.
99 p. : il.
Tese (Doutorado) – Universidade Federal de Lavras, 2008.
Orientador: Júlio Sílvio de Sousa Bueno Filho.
Bibliografia.
1. Algoritmo Metropolis-Hastings. 2. MCMC. 3. Modelos lineares
generalizados. 4. Modelos mistos. 5. Modelo de limiar. I. Universidade
Federal de Lavras. II. Título.
CDD – 519.542
JOSÉ WALDEMAR DA SILVA
ALGORITMOS PARA MODELOS DE LIMIAR USANDO AS DISTRIBUIÇÕES
ACUMULADAS NORMAL E t DE STUDENT
Tese apresentada à Universidade Federal de Lavras
como parte das exigências do Programa de Pós-
graduação em Estatística e Experimentação Agro-
pecuária, para a obtenção do título de “Doutor”.
APROVADA em 07 de março de 2008
Prof
a
. Dr
a
. Roseli Aparecida Leandro ESALQ/USP
Prof
a
. Dr
a
. Joelma Pereira UFLA
Prof
a
. Dr
a
. Thelma Sáfadi UFLA
Prof. Dr. Daniel Furtado Ferreira UFLA
Prof. Dr. Júlio Sílvio de Sousa Bueno Filho
UFLA
(Orientador)
LAVRAS
MINAS GERAIS-BRASIL
Aos meus pais, Gregório (in memoriam) e Maria Edith,
por serem exemplos de determinação, sempre acreditarem
na educação de seus filhos e incentivá-los.
E aos meus irmãos, Valdete, José Valdo,
Valter, Valdeir e Valdinei,
pelo carinho e
amizade,
Ofereço.
À minha esposa Elizângela, pela amizade, compreensão e amor,
Dedico.
AGRADECIMENTOS
A Deus, por tudo.
À Universidade Federal de Lavras (UFLA) e ao Programa de Pós-Graduação em
Estatística e Experimentação Agropecuária por propiciar a realização deste trabalho.
Ao professor Júlio Sílvio de Sousa Bueno Filho pela orientação e ensinamentos.
Aos professores membros da banca examinadora, pelas críticas, sugestões e co-
laboração.
Aos professores do Departamento de Ciências Exatas (DEX), pelos ensinamentos
fundamentais para a minha formação.
A todos os funcionários do DEX, pela atenção e eficiência com que nos atende-
ram.
Ao colega de doutorado, Janser Moura Pereira, pela amizade e companheirismo
desde o início da caminhada no curso de graduação.
Aos colegas de curso pós-graduação, pelo incetivo, companheirismo e frater-
nidade, em especial aos amigos Imaculada, Claudiney, Verônica, Andréia, Denismar,
Fabyano, Cirilo, Nádia, Rosi e Renato.
À Renata, pela amizade e por se dispor a me acolher sempre.
Aos amigos Célio e Célia, por me acolherem sempre prontamente em sua casa.
Aos amigos Airton e Andréia por tudo e pelos momentos agradáveis que vivenciamos
em Teresina-PI. Aos amigos, Isaías, Luciana e Larissa, Flávio e Carminha por me propi-
ciarem momentos de descontração e alegria.
A todos os meus familiares, que contribuíram para a minha formação humana,
acreditaram e colaboraram para o meu sucesso profissional e pessoal.
A todos aqueles que, direta ou indiretamente, contribuíram para a conclusão desta
etapa, muito importante em minha vida, o meu sincero agradecimento.
SUMÁRIO
LISTA DE TABELAS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . i
LISTA DE FIGURAS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii
RESUMO . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv
ABSTRACT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi
1 INTRODUÇÃO . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
2 REFERENCIAL TEÓRICO . . . . . . . . . . . . . . . . . . . . . . . . 4
2.1 Algoritmo de Sörensen et al. (1995) - (AC-normal) . . . . . . . . . . . 8
2.2 Algoritmo de Cowles (1996) e Kizilkaya et al. (2003) - (MC-normal) . 13
2.3 Algoritmo de Kizilkaya et al. (2003) - (MC-t) . . . . . . . . . . . . . . 15
2.4 Comparação de Modelos . . . . . . . . . . . . . . . . . . . . . . . . . 22
3 MATERIAL E MÉTODOS . . . . . . . . . . . . . . . . . . . . . . . . . 26
3.1 Material . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
3.2 Métodos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
4 RESULTADOS E DISCUSSÃO . . . . . . . . . . . . . . . . . . . . . . 34
4.1 Resultados metodológicos . . . . . . . . . . . . . . . . . . . . . . . . 34
4.1.1 Algoritmo adaptado de AC (AC-t) . . . . . . . . . . . . . . . . . . . 34
4.1.2 Algoritmo adaptado de NC - (NC-normal) . . . . . . . . . . . . . . . 35
4.1.3 Algoritmo adaptado de NC - (NC-t) . . . . . . . . . . . . . . . . . . 40
4.2 Resultados do ajuste para os três algoritmos . . . . . . . . . . . . . . . 44
4.2.1 Algoritmo AC (t e normal) . . . . . . . . . . . . . . . . . . . . . . . 44
4.2.2 Algoritmo MC (t e normal) . . . . . . . . . . . . . . . . . . . . . . . 50
4.2.3 Algoritmo NC (t e normal) . . . . . . . . . . . . . . . . . . . . . . . 54
4.2.4 Análise dos contrastes de interesse . . . . . . . . . . . . . . . . . . . 61
4.3 Comparações com a ANAVA e considerações gerais . . . . . . . . . . 62
5 CONCLUSÕES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
REFERÊNCIAS BIBLIOGRÁFICAS . . . . . . . . . . . . . . . . . . . . 67
ANEXO . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
LISTA DE TABELAS
1 Escala hedônica de nove pontos. . . . . . . . . . . . . . . . . . . 4
2 Formas de análise encontradas na literatura para o problema de
dados categorizados em escalas subjetivas. . . . . . . . . . . . . . 21
3 Escala para interpretação do fator de Bayes. . . . . . . . . . . . . 24
4 Notas atribuídas por 36 provadores à cor da banana da terra desi-
dratada utilizando três concentrações de sacarose. . . . . . . . . . 27
5 Esquema de reclassificação dos dados em cinco categorias. . . . . 27
6 Valores de Burn-in (B) e Thinning (T) obtidos por meio do diag-
nóstico de Raftery & Lewis, para as médias de tratamentos (β),
variância de blocos (σ
2
u
), variância residuals (δ
2
) e graus de liber-
dade (ν), nos algoritmos AC, MC e NC segundo a distribuição da
variável latente. . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
7 Autocorrelação entre as estimativas para β
1
e para γ
3
no caso em
que a função de ligação foi a distribuição normal acumulada e para
β
2
e γ
2
quando a função de ligação foi a distribuição t de Student
acumulada, no algoritmo AC. . . . . . . . . . . . . . . . . . . . . 49
8 Autocorrelação entre as estimativas de β
1
e de γ
4
no caso em que
a função de ligação foi a distribuição normal acumulada e de β
1
e
γ
3
quando foi usada a distribuição t de Student acumulada como
função de ligação, no algoritmo MC. . . . . . . . . . . . . . . . . 53
9 Autocorrelação entre as estimativas de β
1
e de γ
2
no caso em que a
função de ligação foi a distribuição normal acumulada e de β
1
e de
γ
3
quando foi usada a distribuição t de Student acumulada como
função de ligação, no algoritmo NC. . . . . . . . . . . . . . . . . 58
i
10 HPDs para funções dos parâmetros associadas à regressão linear
(C
1
β), ao desvio de regressão (C
2
β) e à correlação intraclasse (ρ). 62
11 Quadro-resumo da análise da variância para a variável cor sem
transformação, com transformação logarítmica e com transforma-
ção Box-Cox. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
ii
LISTA DE FIGURAS
1 Funções de distribuição normal, t de Student e logística de acordo
com o quantil. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2 Cadeias para médias de tratamentos no algoritmo AC e nas duas
funções de ligação normal de t de Student, acumuladas. . . . . . . 45
3 Cadeias para os parâmetros de limiar no algoritmo AC e modelos
normal e t de Student. . . . . . . . . . . . . . . . . . . . . . . . . 46
4 Cadeia e distribuição das estimativas de ν no algoritmo AC. . . . . 49
5 Cadeias para médias de tratamentos no algoritmo MC e nas duas
funções de ligação normal de t acumuladas. . . . . . . . . . . . . 51
6 Cadeias para os parâmetros de limiar no algoritmo MC e modelos
normal e t de Student. . . . . . . . . . . . . . . . . . . . . . . . . 52
7 Distribuição das estimativas de ν no algoritmo MC. . . . . . . . . 54
8 Cadeias para médias de tratamentos no algoritmo NC e nas duas
funções de ligação normal de t de Student, acumuladas. . . . . . . 56
9 Cadeias para os parâmetros de limiar no algoritmo NC e modelos
normal e t. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
10 Cadeia e histograma para ν no algoritmo NC. . . . . . . . . . . . 59
11 Cadeias para os contrastes associados à regressão linear e ao des-
vio de regressão no algoritmo AC e funções de ligação normal e t
de Student, acumuladas. . . . . . . . . . . . . . . . . . . . . . . 63
12 Gráfico de resíduos a partir do modelo ajustado para a variável
nota atribuída à cor da banana da terra desidratada, sem transfor-
mação (a), com transformação logarítmica (b) e com transforma-
ção de Box-Cox (c). . . . . . . . . . . . . . . . . . . . . . . . . . 65
iii
RESUMO
Silva, José Waldemar da. Algoritmos para modelos de limiar usando as dis-
tribuições acumuladas Normal e t de Student. Lavras: UFLA, 2008. 99 p.
Tese (Doutorado em Estatística e Experimentação Agropecuária) - Universidade
Federal de Lavras, Lavras, MG.
*
Dados categorizados podem ser modelados por meio de uma variável la-
tente (L) com escala contínua que relaciona a nota de uma categoria à pertinên-
cia em intervalo correspondente. Modela-se então a esperança de L por meio de
uma função de ligação que é uma diferença entre valores de uma função acumu-
lada nos limites de intervalo. Isto permite o estabelecimento de modelos lineares
para a variável latente e propicia maior simplicidade computacional e facilidade
de interpretação. O amostrador de Gibbs (GS) foi o primeiro algoritmo proposto
(Albert & Chib, 1993; algoritmo AC) para a obtenção de amostras das distribui-
ções conjuntas a posteriori, usando a distribuição t de Student, para o modelo
fixo. Foi constatada, no entanto, forte autocorrelação nas cadeias dos parâmetros
do modelo linear e nas dos parâmetros de limiar (limites dos intervalos para L).
Visando minimizar este problema surgiram dois tipos de propostas de amostra-
gem conjunta da variável latente e dos parâmetros de limiar usando algoritmos do
tipo Metropolis-Hastings. A primeira usa como geradora de candidatos a distri-
buição normal truncada e está aqui implementada em um algoritmo modificado de
Cowles (1996) (algoritmo MC). A segunda usa a distribuição Dirichlet e está no
algoritmo modificado de Nandran & Chen (1996) (algoritmo NC). Neste trabalho
as três metodologias para a modelagem da variável latente foram adaptadas para
a análise de modelos mistos gerais. Os três algoritmos foram estudados quanto à
rapidez de convergência em um experimento real de análise sensorial. Neste ex-
perimento, foi estudado o efeito de três concentrações de sacarose no processo de
desidratação da banana da terra. Foram atribuídas notas para a cor do produto em
uma escala hedônica de nove pontos. Quanto à convergência, a superioridade do
algoritmo NC foi evidente para as duas funções de ligação utilizadas. O modelo
com distribuição t de Student foi o melhor para os algoritmos AC e MC. Para o
caso do algoritmo NC o resultado não discrimina os modelos, o que deve favorecer
o modelo normal. Os resultados são semelhantes aos registrados por análises da
variância usuais, cujas pressuposições são, no entanto, violadas. Os algoritmos es-
tão implementados em R e permitem analisar quaisquer estruturas de covariância
para os efeitos aleatórios.
Palavras-chave: Algoritmo Metropolis-Hastings;MCMC; modelos lineares gene-
ralizados; modelos mistos; modelos de limiar.
*
Orientador:Júlio Sílvio de Sousa Bueno Filho - UFLA.
iv
ABSTRACT
Silva, José Waldemar da. Algorithms for threshold models using Normal and
Student-t cumulative distributions. Lavras: UFLA, 2008. 99 p. Thesis (Doc-
torate in Statistics and Agricultural Experimentation) - Universidade Federal de
Lavras, Minas Gerais, Brazil.
*
Discrete ordinal data can be modeled by a latent variable (L) whose do-
main is a set of real intervals with bijecting relation between cathegory and be-
longing to the respective interval. The expectation of L is then modeled using a
link function that is a difference of cumulative distributions. This methodology
has the advantage of establishing linear models for the latent variable with great
computing simplification and interpretability. Gibbs Sampling (GS) was the first
proposed algorithm to sample from the posterior in these situations (Albert & Chib,
1993; algorithm AC). This algorithm allows Student’s t distribution, for fixed mo-
dels. However, strong autocorrlations where found in Marcov chains from both
of the linear predictors and threshold parameters. Two kinds of proposals, both
based on joint sampling of L and threshold parameters via Metropolis-Hastings
type algorithms, are described to minimize this problem. The first uses a trun-
cated normal distribution as a candidate generating function and is implemented
here in an algorithm modified from Cowles (1996) (MC algorithm). The second
uses a Dirichlet distribution and is implemented in the algorithm modified from
Nandran & Chen (1996) (NC algorithm). In this work, the three methodologies
were adapted to a general mixed model analysis. Convergence readiness in a real
experiment was compared for the three algorithms. A sensorial analysis was used
as an example. In this experiment, the effect of three concentrations of sacharosis
in the dehidrating a banana variety were studied. Grades to the product color were
given in a nine point hedonic scale. Regarding convergence, NC was the best algo-
rithm for both normal and t distributions. Student’s t was the more probable model
using SAGK and MC algorithms, but for NC algorithm there was no evidence of
difference among t and normal models (result that favors normal model). Resul-
ting analysis were similar to usual ANOVA approximation, although violating its
assumptions. Algorithms are implemented in R and allow to analyze a completely
general covariance structure for random effects.
Key-words: MCMC; Metropolis-Hastings Algorithm; generalized linear models;
mixed models; threshold models.
*
Adviser: Júlio Sílvio de Sousa Bueno Filho - UFLA.
v
1 INTRODUÇÃO
Muitos atributos na experimentação, são avaliados por meios subjetivos.
Juízes, técnicos, provadores ou consumidores são escolhidos ou designados para
atribuir um conceito ao material experimental de acordo com suas preferências ou
de acordo com algum regulamento.
Em geral, para a obtenção dos dados adota-se um conjunto de categorias
dotadas de certa ordem nas quais os produtos possam ser classificados do pior ao
melhor, de acordo com o julgamento de cada avaliador, diferentemente dos casos
em que se tem uma variável contínua ou na forma de contagem.
Ainda é muito comum o uso desta metodologia sem a preocupação com a
validade dos testes aplicados, apesar de os dados descritos não seguirem as pres-
suposições da análise de variância.
Outro ponto a ser destacado é a dificuldade ou impossibilidade encontrada,
em muitos casos, para quantificar uma característica por meio de medição. Podem
ser citados alguns exemplos típicos para esta situação como a seleção para algu-
mas características morfológicas de touros, no melhoramento genético de gado de
corte e, ou, leiteiro, tais como escore corporal e a classificação do animal quanto ao
padrão racial, a avaliação do grau de dificuldade do parto de matrizes reprodutoras
(vacas) e o grau de satisfação do consumidor com relação ao sabor de um produto
alimentício a ser lançado no mercado. Vários outros exemplos poderiam ser men-
cionados, dada a grande quantidade de características de importância econômica
existentes que são subjetivas, categóricas e ordinais.
A literatura recente sobre análise de dados categorizados subjetivamente,
sugere o ajuste de modelos com parâmetros threshold ou de limiar, no qual, uma
variável latente, com distribuição contínua, é especificada de forma que uma res-
1
posta é observada em uma dada categoria se o valor desta variável está entre os
limites que definem tal categoria.
A variável latente facilita a implementação computacional, pois, pode ser
modelada de forma linear. Esta variável pode ser interpretada como o estímulo
que deve ser aplicado a um indivíduo, objeto ou produto, para que este possa ser
classificado em uma das K categorias. A área da análise sensorial dos alimentos é
um exemplo disto pois, muitas características são difíceis ou impossíveis de serem
quantificadas por meio de medida objetiva. Nestes casos, cada provador ou juiz
classifica um produto em alguma categoria, por meio do seu sabor, cor, textura,
etc. Estas variáveis são dotadas de uma escala contínua, porém, a execução do
experimento torna-se mais fácil quando são fragmentadas em categorias.
Wright (1934) foi o primeiro a propor esta metodologia, no estudo da he-
rança do número de dedos de Cavia porcellus (porquinhos da Índia). Albert &
Chib (1993) propuseram a utilização da amostragem Gibbs para obter aproxima-
ção das distribuições marginais dos parâmetros. Sörensen et al. (1995) utilizaram
a mesma metodologia, porém na área da Genética Quantitativa Animal em que o
Modelo Linear Generalizado (MLG) correspondente usava uma função de ligação
normal acumulada.
A forte autocorrelação entre as estimativas, principalmente dos parâmetros
de limiar, obtidos por meio da amostragem Gibbs, usando esta metodologia, é
evidente e foi destacada por Sörensen et al. (1995) e Sörensen & Gianola (2004).
Este fato obriga a obtenção de uma cadeia muito grande para que se tenha uma boa
representação da distribuição marginal do parâmetro. Propostas de algoritmos para
contornar ou amenizar este problema são dados por Cowles (1996) e Nandram &
Chen (1996).
A escolha de um bom modelo e a conseqüente maior qualidade no processo
2
de inferência podem ser quantificados por alguns indicadores como a maior razão
de verossimilhanças marginais dos modelos que se pretende comparar, a maior
acurácia nas estimativas ou a melhor interpretabilidade de seus parâmetros.
O algoritmo para a análise de modelos mistos com função de ligação
dada por uma distribuição normal acumulada foi desenvolvida por Sörensen et
al. (1995) e é apresentada na seção 2.1. O mesmo algoritmo porém, para modelos
fixos com função de ligação t de Student acumulada, é apresentado por Albert &
Chib (1993). Contudo, outro algoritmo utilizando a amostragem conjunta dos pa-
râmetros threshold e da variável latente foi proposto por Cowles (1996) e adaptado
por Kizilkaya et al. (2003) para modelos mistos no melhoramento animal. Este
algoritmo é apresentado nas seções 2.2, com a distribuição normal padrão como
funçao de ligação, e 2.3, com t de Student acumulada. Ainda outra melhoria nes-
tes algortimos é a proposta de Nandran & Chen (1996) que visa reparametrizar os
threshold e a variável latente.
Este trabalho teve como objetivo adaptar o algoritmo de Albert & Chib
(1993) e de Nandran & Chen (1996) para os modelos mistos, além de usar a distri-
buição t de Student como função de ligação no segundo e contrastar a velocidade
de convergência dos algoritmos propostos por Albert & Chib (1993) - AC, Co-
wles (1996) - MC e Nandram & Chen (1996) - NC e adaptados. A aplicação da
metodologia será feita em um experimento de análise sensorial de alimentos.
Um objetivo adicional foi o de implementar estas metodologias em rotinas
flexíveis no pacote estatístico R, que permitam análise de modelos mistos em si-
tuações frequentes e com generalidade suficiente para incluir estruturas de efeitos
aleatórios correlacionados. Tais rotinas estão disponibilizadas em anexo.
3
TABELA 1: Escala hedônica de nove pontos.
Categorias E. H. de 9 Pontos
Desgostei muitíssimo 1
Desgostei muito 2
Desgostei 3
Desgostei moderadamente 4
Nem gostei nem desgostei 5
Gostei moderadamente 6
Gostei 7
Gostei muito 8
Gostei muitíssimo 9
Fonte: Dutcosky (1996).
2 REFERENCIAL TEÓRICO
Alguns atributos na análise sensorial de alimentos são avaliados por meio
da escala hedônica de nove pontos, conforme Tabela 1.
Porém, a existência de um grande número de categorias para a classificação
pode dificultar a decisão do provador quanto a qual nota atribuir ao produto. Além
disso, o acréscimo do número de parâmetros de limiar estimados pode dificultar a
convergência do algoritmo e sob o aspecto da representatividade, Kizilkaya et al.
(2003) sugerem o agrupamento de categorias com poucas observações.
O conhecimento prévio sobre um fenômeno e a informação fornecida por
experimentos podem ser combinados por meio do Teorema de Bayes, usando a
seguinte expressão:
p(θ|Y ) ∝ p(θ)p(Y |θ),
em que θ é o vetor de parâmetros, Y o vetor de observações, p(θ|Y ) a distribuição
conjunta a posteriori e p(θ) a distribuição a priori.
As inferências são realizadas a partir da distribuição conjunta a posteriori
4
que em geral tem forma complexa e se torna difícil a avaliação das quantidades
de interesse. O problema a ser resolvido é o cálculo da esperança a posteriori de
alguma função g(θ),
E[g(θ)|Y ] =
g(θ)p(θ|Y )dθ, (1)
A partir de uma amostra θ
1
, θ
2
, ··· θ
n
da distribuição a posteriori p(θ|Y )
a integral em (1) é aproximada por meio da média,
E[g(θ)|Y ] =
1
n
n
i=1
g(θ
i
).
Geralmente este problema é transferido para o caso univariado e assim,
ele seria solucionado calculando-se a distribuição marginal para cada parâmetro.
No entanto, a abordagem univariada não contorna a dificuldade de integração para
a obtenção da esperança a posteriori. Métodos numéricos como a amostragem
Gibbs podem ser usados para aproximar a integral em (1) quando as distribui-
ções condicionais completas a posteriori possuem formas algébricas fechadas (das
quais se pode tomar amostras). Caso isso não ocorra, deve-se lançar mão de ou-
tras formas de amostragem, em geral baseadas no algoritmo Metropolis-Hastings
(Metropolis et al., 1953 e Hastings, 1970).
Os modelos lineares generalizados (GLM) podem ser usados para estabe-
lecer uma relação entre uma variável resposta categorizada ordinal e variáveis pre-
ditoras, sendo que, a função ligadora é alguma função de distribuição. Supondo
que existam n variáveis, classificadas em K categorias ordenadas Y
1
, Y
2
, ..., Y
n
,
em que Y
i
é tal que p
ik
= P (Y
i
= k) para k = 1, 2, . . . , K e
K
k=1
p
ik
= 1, i =
1, . . . , n e assumindo que as probabilidades acumuladas possam ser modelada por
5
um GLM, têm-se:
P (Y
i
≤ k) = G(γ
k
− w
i
θ), k = 1, 2, ..., K (2)
em que w
i
é o vetor linha de incidência, ligando θ à i-ésima observação, G é uma
função de distribuição de uma variável contínua que pode assumir várias formas,
θ é o vetor de efeitos fixos e aleatórios e γ é o vetor de parâmetros threshold que
divide a reta real em K intervalos disjuntos, da seguinte forma:
(γ
0
, γ
1
); [γ
1
, γ
2
); ··· ; [γ
K−1
, γ
K
);
em que γ
0
= −∞ e γ
K
= +∞
Em geral adota-se uma origem para a variável latente fixando um dos pa-
râmetros de limiar ou threshold e em geral, escolhe-se γ
∗
1
= 0. A motivação para
γ
∗
1
= 0 é dada pela tranformação γ
∗
k
= γ
k
− γ
1
. Assim, cada γ
∗
k
será a distância
entre γ
k
e γ
1
. A variável latente transformada será obtida por meio da expressão
L
∗
i
= L
i
− γ
1
, i = 1, ..., n e assim,
y
i
= k se γ
k−1
− γ
1
≤ L
i
− γ
1
≤ γ
k
− γ
1
Os novos parâmetros (θ
∗
i
), de efeitos fixos e aleatórios, serão θ
∗
i
= θ
i
− γ
1
, i =
0, 1, ..., m.
Na escala da variável latente, L
∗
= W θ
∗
+ em que, W é a matriz de
incidência, θ
∗
o vetor de efeitos e o vetor de resíduos. Assumindo distribuição
normal para a variável latente (L
∗
), então
(L
∗
|θ
∗
, σ
2
e
) ∼ N(W θ
∗
, Iσ
2
e
).
Sem perda de generalidade a variância residual pode ser constante, adotando-se
6
por exemplo, σ
2
e
= 1 (Sörensen & Gianola, 2004).
A distribuição conjunta a posteriori, de todos os parâmetros e da variável
latente resulta em:
p(θ
∗
, L
∗
, γ
∗
|Y ) ∝ p(Y |θ
∗
, L
∗
, γ
∗
)p(θ
∗
, γ
∗
) = p(Y |L
∗
, γ
∗
)p(θ
∗
, L
∗
, γ
∗
)
= p(Y |L
∗
, γ
∗
)p(θ
∗
, L
∗
, γ
∗
) = p(Y |L
∗
, γ
∗
)p(L|θ
∗
)p(θ
∗
, γ
∗
)
= p(Y |L
∗
, γ
∗
)
n
i=1
p(L
∗
i
|θ
∗
)
p(θ
∗
, γ
∗
) (3)
A distribuição p(Y |L
∗
, γ
∗
) em (3) é degenerada pois, a probabilidade de
uma determinada observação pertencer a uma determinada categoria, dado o valor
da variável latente e os limiares, é completamente especificada. Conforme notação
de Albert & Chib (1993), p(Y |L
∗
, γ
∗
) pode ser escrita como
p(Y |L
∗
, γ
∗
) =
n
i=1
K
k=1
I
[γ
∗
k−1
,γ
∗
k
)
(L
∗
i
)I(Y
i
= k)
.
Albert & Chib (1993) apresentam um algoritmo para a modelagem de da-
dos categorizados, no qual a função de ligação (G) em (2) corresponde à distri-
buição t de Student acumulada. O mesmo algortimo é implementado por Sörensen
et al. (1995) em que a função de ligação em (2) é a distribuição normal padrão
acumulada considerando o modelo misto (efeitos fixos e aleatórios no modelo
para a variável latente). Nesta modelagem, um parâmetro de limiar ou threshold
(γ
∗
k
) tem distribuição condicional completa uniforme com domínio no intervalo
[max(L
∗
|Y = k); min(L
∗
|Y = k + 1)] em que L
∗
é o vetor de variáveis latentes.
O algoritmo desenvolvido por Sörensen et al. (1995) é descrito a seguir.
7
2.1 Algoritmo de Sörensen et al. (1995) - (AC-normal)
Supondo que a função G em (2) é dada por uma distribuição normal padrão
acumuada (Φ(·)), então o modelo amostral é
p(Y
i
= k|θ
∗
, γ
∗
) = Φ
γ
∗
k
− w
i
θ
∗
− Φ
γ
∗
k−1
− w
i
θ
∗
. (4)
e usando uma variável latente L
∗
podemos reescrever (2) dados os parâmetros de
limiar γ
∗
, θ
∗
e w
i
como:
Y
i
= k se γ
∗
k−1
≤ L
i
≤ γ
∗
k
L
∗
i
∼ N(w
i
θ
∗
, 1) (5)
com i = 1, 2,..., n e k = 1, 2,..., K.
Para a variável latente L
∗
, tem-se o modelo linear L
∗
= W θ
∗
+ , em
que W = [X|Z]. com X e Z correspondendo às matrizes de delineamento para
os efeitos fixos e aleatórios, respectivamente; é o vetor de erros aleatórios e
θ
∗
= (β
, u
) é o vetor de efeitos fixos (β) e aleatórios (u). Quando não se
tem efeitos aleatórios o modelo terá sua estrutura simplificada, com W = X e
θ
∗
= β
.
Assumindo que os elementos de Y são condicionalmente independentes
dado θ
∗
, a densidade conjunta a posteriori de todos os parâmetros, incluindo a
variável latente L
∗
, é dada como segue:
p(θ
∗
, γ
∗
, L
∗
|Y ) ∝ p(Y, L
∗
, θ
∗
, γ
∗
) = p(Y, L
∗
|θ
∗
, γ
∗
)p(θ
∗
, γ
∗
)
= p(Y |L
∗
, θ
∗
, γ
∗
)p(L
∗
|θ
∗
, γ
∗
)p(θ
∗
, γ
∗
) (6)
∝
n
i=1
exp(−
1
2
(L
∗
i
− w
i
θ
∗
)
2
)
I
[γ
∗
k−1
,γ
∗
k
)
(L
∗
i
)φ(θ
∗
; 0, V )
8
A informação a priori sobre (θ
∗
, γ) pode ser dada por:
p(θ
∗
, γ) ∝ p(θ
∗
) ∝ exp
−
τ
2
θ
∗
θ
∗
a qual será vaga se o parâmetro de precisão τ assumir um valor pequeno, por
exemplo, τ = 0, 001 ou em (6) fazendo V = 1000I
m
, em que I
m
é uma matriz
identidade m × m.
Quando se tem um modelo linear misto, com θ
= (β
, u
), u é especifi-
cado por
u|σ
2
u
∼ N
∅, A
q
σ
2
u
em que A é a matriz de correlação entre os efeitos aleatórios. Em modelos de
Genética Quantitativa, por exemplo, A é um múltiplo da matriz de parentesco
(Kempthorne, 1966).
Assim, a matriz V pode ser dada por
V =
1000I
p
∅
∅ σ
2
u
A
q
em que σ
2
u
é a variância dos efeitos aleatórios, p q correspondem, respectivamente,
à quantidade de efeitos fixos e à quantidade de efeitos aleatórios, I
p
é uma matrizes
identidade p × p e I
q
é a matriz de correlação entre os elementos de u, q × q.
No modelo misto,
p(θ
∗
, γ) ∝ p(θ
∗
) ∝
1
σ
2
u
(
q
2
)
exp
−
1
2
θ
∗
V
−1
θ
∗
.
9
A priori para σ
2
u
pode ser dada por uma Gama Inversa,
p
σ
2
u
=
σ
2
u
−(a+1)
exp
−b
σ
2
u
, (7)
em que a e b são os parâmetros desta distribuição a priori, a qual tem desejáveis
propriedades de conjugação com os modelos normais.
E assim, a distribuição conjunta em (6) é reescrita como segue:
p(θ
∗
, γ
∗
, L
∗
, σ
2
u
|Y ) ∝ p(Y, L
∗
, θ
∗
, γ
∗
) = p(Y, L
∗
|θ
∗
, γ
∗
)p(θ
∗
, γ
∗
)
= p(Y |L
∗
, θ
∗
, γ
∗
)p(L
∗
|θ
∗
, γ
∗
)p(θ
∗
, γ
∗
)
∝
n
i=1
exp(−
1
2
(L
∗
i
− w
i
θ
∗
)
2
)
I
[γ
∗
k−1
,γ
∗
k
)
(L
∗
i
)
× φ(θ
∗
; 0, V )p(σ
2
u
) (8)
As inferências sobre cada parâmetro podem ser feitas com base em amos-
tras da distribuição conjunta a posteriori (Gelfand & Smith, 1990), calculadas a
partir de (6), se o interesse for no modelo fixo ou a partir de (8) se o interesse
for no modelo misto, usando os métodos Monte Carlo via Cadeias de Markov
(MCMC).
A metodologia de amostragem das distribuições condicionais foi original-
mente desenvolvida por Sörensen et al. (1995) a partir de (8) e está apresentada
também, de forma bem didática, em Sörensen & Gianola (2004). Nos parágrafos
que se seguem estamos reproduzindo a maior parte desta apresentação.
É preciso notar que a distribuição p(Y |L
∗
, γ
∗
) na expressão em (8) é de-
generada, pois o conhecimento de L
∗
leva ao conhecimento exato de Y (Albert &
10
Chib, 1993) e sua expressão é dada por:
p (Y |L
∗
, γ
∗
) =
n
i=1
K
k=1
I
[
γ
∗
k−1
,γ
∗
k
)
(L
∗
i
)
I (Y
i
= k)
(9)
A partir da expressão anterior pode-se observar que a distribuição condi-
cional completa de γ
∗
k
dado γ
∗
−k
, θ
∗
e L
∗
, em que γ
∗
−k
representa todos os parâme-
tros de limiar exceto γ
k
, é dada por:
P
γ
∗
k
|γ
∗
−k
, L
∗
, Y
∝ (10)
∝
n
i=1
I (Y
i
= k) I
[
γ
∗
k−1
,γ
∗
k
)
(L
∗
i
) + I (y
i
= k + 1) I
[
γ
∗
k
,γ
∗
k+1
)
(L
∗
i
)
a qual pode ser vista como uma distribuição uniforme, conforme expressão a seguir
P
γ
∗
k
|γ
∗
−k
, L
∗
, Y
=
1
a − b
, a > b (11)
em que a = min
(L
∗
|Y = k + 1) , γ
∗
k+1
e b = max
(L
∗
|Y = k) , γ
∗
k−1
.
Com relação à variável latente, observa-se que L
∗
i
tem distribuição condi-
cional completa a posteriori dada por:
P
L
∗
i
|θ
∗
, Y
i
= k, γ
∗
k−1
, γ
∗
k
=
φ
L
∗
i
w
i
θ
∗
, 1
Φ
γ
∗
k
− w
i
θ
∗
− Φ
γ
∗
k−1
− w
i
θ
∗
,
(12)
com γ
∗
k−1
< L
∗
i
≤ γ
∗
k
, Φ(.) representando a distribuição normal padrão acumulada
e φ
L
∗
i
(w
i
θ
∗
, 1) a densidade normal com média w
i
β
∗
e variância unitária.
O vetor de parâmetros θ
∗
tem distribuição condicional completa a posteri-
11
ori dada por uma normal multivariada, ou seja,
θ
∗
|Y, L
∗
∼ N
B
−1
W
L
∗
, B
−1
(13)
em que B = V
−1
+ W
W .
Para a variância dos efeitos aleatórios (σ
2
u
), a distribuição condicional
completa será
p
σ
2
u
|θ
∗
, L
∗
, Y
=
σ
2
u
−(
q
2
+a+1)
exp
−1
2σ
2
u
u
u + 2b
, (14)
que é o núcleo de uma distribuição Gamma Inversa com parâmetros como a seguir:
σ
2
u
|θ
∗
, L
∗
, Y ∼ GI
q + 2a
2
,
u
u + 2b
2
(15)
A amostragem Gibbs, segundo Sörensen et al. (1995), pode ser implemen-
tada a partir das distribuições condicionais completas dadas em (13), (15), (12) e
(11), não necessariamente nesta ordem.
Um problema apontado por alguns autores como Sörensen & Gianola (2004),
Cowles (1996), Nandran & Chen (1996) e Kizilkaya et al. (2003) é que nos algo-
ritmos apresentados por Albert & Chib (1993) e por Sörensen et al. (1995) existe
uma forte autocorrelação entre os valores obtidos por meio da amostragem Gibbs,
principalmente para os parâmetros de limiar. Segundo Nandran & Chen, (1996)
dificuldade ou lentidão na convergência para esses parâmetros acarreta, por con-
seqüência, dificuldade de convergência para θ.
Algumas propostas de algoritmos para contornar ou amenizar este pro-
blema são dados por Cowles (1996), Kizilkaya et al. (2003) e Nandram & Chen
(1996).
12
Para contornar ou amenizar problemas de convergência Cowles (1996)
apresenta uma análise em que amostras da distribuição conjunta a posteriori são
obtidas por meio da amostragem Gibbs com um “passo” Metropolis-Hastings para
obter amostras dos parâmetros threshold e da variável latente a partir da forma fa-
torada p(γ
∗
, L
∗
|θ
∗
, Y ) = p(L
∗
|γ
∗
, θ
∗
, Y )p(γ
∗
|θ
∗
, Y ). Nesta análise a função de
ligação utilizada em (2) foi a distribuição normal padrão acumulada e um modelo
fixo para a variável latente. De acordo com o método da composição, apresentado
por Devroye (1986, citado por Gelfand & Smith, 1990), Chib (2001) e Sörensen &
Gianola (2004), a distribuição de (γ
∗
, L
∗
|θ
∗
, Y ) pode ser escrita como o produto
das distribuições condicionais (γ
∗
|θ
∗
, Y ) e (L
∗
|γ
∗
, θ
∗
, Y ), conforme teorema de
Bayes em que (L
∗
|γ
∗
, θ
∗
, Y ) é obtida a partir da posteriori conjunta e (γ
∗
|θ
∗
, Y )
é obtida da integração de (γ
∗
, L
∗
|θ
∗
, Y ) em relação à L
∗
. Conforme Liu et al.
(1994) agrupando os elementos a serem amostrados, em geral, tem-se aumento na
eficiência do amostrador de Gibbs. A candidata utilizada para gerar valores de
γ
∗
k
nesta metodologia foi a distribuição normal, truncada por outros dois parâme-
tros threshold adjacentes, com média dada pelo threshold amostrado na iteração
anterior e uma variância fixa.
Seguindo a metodologia apresentada por Cowles (1996), Kizilkaya et al.
(2003) apresentam um modelo misto para a variável latente associada à variável
escore para facilidade de parto, utilizando como funções de ligação as distribuições
normal e t de Student acumuladas.
2.2 Algoritmo de Cowles (1996) e Kizilkaya et al. (2003) - (MC-normal)
Visando melhorar o processo de amostragem, é apresentada a seguir, a
proposta de Kizilkaya et al. (2003), adaptada de Cowles (1996). Segundo esta
metodologia, a atualização dos parâmetros de limiar não é feita de forma isolada,
13
como proposto por Albert & Chib (1993), mas de forma a aceitar todo o vetor
formado por estes parâmetros.
A distribuição para (γ
∗
, L
∗
|θ
∗
, Y ) pode ser fatorada como o produto das
distribuições condicionais (γ
∗
|θ
∗
, Y ) e (L
∗
|γ
∗
, θ
∗
, Y ). A distribuição condicional
para (γ
∗
|θ
∗
, Y ) é obtida a partir da integral de (γ
∗
, L
∗
|θ
∗
, Y ) em relação a L
∗
:
p (γ
∗
|θ
∗
, Y ) ∝
Y
i
=2
Φ
γ
∗
2
− w
i
θ
∗
− Φ
−w
i
θ
∗

×
Y
i
=3
Φ
γ
∗
3
− w
i
θ
∗
− Φ
γ
∗
2
− w
i
θ
∗

···
··· ×
Y
i
=K−1
Φ
γ
∗
K−1
− w
i
θ
∗
− Φ
γ
∗
K−2
− w
i
θ
∗

×
Y
i
=K
1 − Φ
γ
∗
K−1
− w
i
θ
∗

(16)
e a distribuição de (L
∗
|γ
∗
, θ
∗
, Y ) resulta em:
L
∗
|γ
∗
, θ
∗
, Y
i
= k ∼ N
w
i
θ
∗
, 1
, (17)
com γ
∗
k−1
< L
∗
i
≤ γ
∗
k
.
Observe que para a amostragem de γ
∗
é necessário o uso do algoritmo
Metropolis-Hastings pois estes parâmetros não têm uma distribuição distribuição
da qual se possa amostrar diretamente.
Cowles (1996) propôs usar a distribuição normal, para a geração de can-
didatos para um determinado threshold, truncada por outros dois adjacentes. Isto
é,
γ
∗
k
|γ
∗
k−1,j
, γ
∗(j−1)
k+1
∼ N
γ
∗(j−1)
k
, σ
2
γ
. (18)
em que j representa a j-ésima iteração e γ
∗
k−1,j
≤ γ
∗
k
< γ
∗(j−1)
k+1
. A probabilidade
14
de aceitação do novo vetor de parâmetros threshold é o mínimo entre 1 e R, em
que
R =
K−1
k=2
Φ

γ
∗(j−1)
k+1
− γ
∗(j−1)
k
/σ
γ
− Φ

γ
∗
k−1,j
− γ
∗(j−1)
k
/σ
γ
Φ

γ
∗
k+1,j
− γ
∗
k,j
/σ
γ
− Φ

γ
∗(j−1)
k−1
− γ
∗
k,j
/σ
γ
×
n
i=1
Φ
γ
∗
Y
i,j
− w
i
θ
∗
− Φ
γ
∗
Y
i−1,j
− w
i
θ
∗
Φ
γ
∗(j−1)
Y
i
− w
i
θ
∗
− Φ
γ
∗(j−1)
Y
i−1
− w
i
θ
∗
(19)
O amostrador de Gibbs pode ser implementado a partir de (13), (15), (18)
e (17).
Segundo Albert & Chib (1993), inferências podem ser sensível à função
de ligação adotada para G em (2). Alternativamente à distribuição normal, Albert
& Chib (1993), Kizilkaya et al. (2003) atribuíram distribuição t de Student para a
variável latente. A distribuição logística investigada por Albert & Chib, (1993) na
análise de dados dicotomizados é outra alternativa (Figura 1).
2.3 Algoritmo de Kizilkaya et al. (2003) - (MC-t)
Em geral a distribuição t de Student é considerada mais plástica para da-
dos binários, mas no caso dos dados categorizados, é provável que tal distribuição
possa assumir formas mais próximas tanto da distribuição normal quanto da logís-
tica, com diferentes gradações dadas pelo número de graus de liberdade. A distri-
buição t de Student aproxima-se da distribuição normal quando o número de graus
de liberdade ν é grande (ν → ∞) e da distribuição logística quando (ν → 0).
Na Figura (1) estão ilustradas as distribuições normal padrão, t de Student com 3
graus de liberdade e logística.
Para a utilização da distribuição t de Student acumulada em (2) é interes-
15
Quantil
Probabilidade Acumulada
Normal
t de Student
Logística
0.0 0.2 0.4 0.6 0.8 1.0
−5 −4 −3 −2 −1 0 1 2 3 4 5
FIGURA 1: Funções de distribuição normal, t de Student e logística de acordo
com o quantil.
sante escrevê-la em dois estágios (Paulino et al., 2003, Sörensen & Gianola 2004),
isto é, como a mistura de normais (21) com a gamma invertida (22) para a vari-
ância. Tendo em vista que a maioria das distribuições utilizadas na modelagem
são da família exponencial, este procedimento propiciará facilidades algébricas e
melhor tratabilidade analítica para as condicionais completas a posteriori.
Utilizando densidade acumulada de uma distribuição t de Student com ν
graus de liberdade a probabilidade de Y pode ser modelada por meio da expressão
(20).
P (Y
i
= k|θ
∗
, ν, γ
∗
) = F
ν
γ
∗
k
− w
i
θ
∗
− F
ν
γ
∗
k−1
− w
i
θ
∗
, (20)
k = 1, 2, ..., K.
Amsotras da distribuição t de Student para a variável latente podem ser
16
obtidas por um processo hierárquico de amostragem de distribuição normal com
variância modificada. Desta forma, amostra-se inicialmente os parâmetros modi-
ficadores λ a partir da distribuição em (22) e em seguida, utilizando estes valores,
amostra-se L de (21).
Os dois estágios para a especificação da distribuição t de Student são indi-
cados nas expressões (21) e (22).
L
∗
i
|θ
∗
, λ
i
∼ N
w
i
θ
∗
,
1
λ
i
, i = 1, 2, ..., n (21)
λ
i
|ν ∼ Gamma
ν
2
,
ν
2
. (22)
Com esta especificação, o modelo amostral (20) pode ser reescrito da se-
guinte maneira:
P (Y
i
= k|θ
∗∗
, λ
i
, ν, γ
∗∗
) = Φ
γ
∗∗
k
− w
i
θ
1
√
λ
i
− Φ
γ
∗∗
k−1
− x
i
θ
1
√
λ
i
, (23)
k = 1, 2, ..., K
Admitindo priori vaga para θ
∗
, conforme especificado em 2.1, a posteriori
conjunta de todos os parâmetros e a variável auxiliar ou latente é:
p(θ
∗
, γ
∗
k
, L
∗
, λ, σ
2
u
, ν|Y ) ∝
n
i=1
φ
L
∗
i
, w
i
θ
∗
,
1
λ
i
I
[
γ
∗
k−1
,γ
∗
k
)
(L
∗
i
)
×
n
i=1
λ
(
ν
2
)
−1
i
exp
−
λ
i
2
ν
× φ (θ
∗
, 0, V ) p
σ
2
u
p (ν) (24)
em que λ = {λ
i
}
n
i=1
.
A matriz V , em (24) é dada na seção 2.1, a priori para σ
2
u
, é como em
17
(7) e a priori para ν é p(ν) = 1/(1 + ν)
2
, isto representa pequena probabilidade
de ocorrência de valores altos para os graus de liberdade, já que, para ν grande a
distribuição t de Student aproxima-se da distribuição normal.
Da mesma forma com que foi descrito para a distribuição normal, a amos-
tragem de γ
∗
e L
∗
é feita de forma conjunta. Quando a função de ligação é a
distribuição t de Student acumulada, a distribuição de (γ
∗
|θ
∗
, Y ) resulta em:
p (γ
∗
|θ
∗
, Y ) ∝
Y
i
=2
Φ
γ
∗
2
− w
i
θ
∗
1
√
λ
i
− Φ
−w
i
θ
∗
1
√
λ
i

×
Y
i
=3
Φ
γ
∗
3
− w
i
θ
∗
1
√
λ
i
− Φ
γ
∗
2
− w
i
θ
∗
1
√
λ
i

···
··· ×
Y
i
=K−1
Φ
γ
∗
K−1
1 − w
i
θ
∗
1
√
λ
i
− Φ
γ
∗
K−2
− w
i
θ
∗
1
√
λ
i

×
Y
i
=K
1 − Φ
γ
∗
K−1
− w
i
θ
∗
1
√
λ
i

e L
∗
tem distribuição normal com a variância modificada por λ como segue:
L
∗
|γ
∗
, θ
∗
, Y
i
= k ∼ N
W θ
∗
,
1
√
λ
i
. (25)
A distribuição candidata para gerar os novos parâmetros de limiar também
pode ser dada por uma normal como em (18) e a probabilidade de aceitação será o
mínimo entre 1 e R, em que R é dado por:
R =
K−1
k=2
Φ

γ
∗(j−1)
k+1
− γ
∗(j−1)
k
/σ
γ
− Φ

γ
∗
k−1,j
− γ
∗(j−1)
k
/σ
γ
Φ

γ
∗
k+1,j
− γ
∗
k,j
/σ
γ
− Φ

γ
∗(j−1)
k−1
− γ
∗
k,j
/σ
γ
×
n
i=1
Φ
γ
∗
Y
i,j
−w
i
θ
∗
1
√
λ
i
− Φ
γ
∗
Y
i−1,j
−w
i
θ
∗
1
√
λ
i
Φ
γ
∗(j−1)
Y
i
−w
i
θ
∗
1
√
λ
i
− Φ
γ
∗(j−1)
Y
i−1
−w
i
θ
∗
1
√
λ
i
. (26)
18
A seguir são apresentadas as distribuições condicionais completas a partir
de (24) para cada parâmetro.
Para o vetor θ a distribuição condicional resulta em:
θ
∗
|Y, L
∗
, λ, γ
∗
∼ N
M
−1
W
R
−1
L
∗
, M
−1
(27)
em que M = W
R
−1
W + V
−1
e R é uma matriz diagonal com seus elementos
dados por λ
−1
i
.
A distribuição condicional completa para λ
i
é dada conforme expressão a
seguir:
p(λ
i
|λ
−i
, L
∗
, θ
∗
, ν) ∝ λ
(
ν+1
2
)
−1
i
exp
−
λ
i
2
L
∗
i
− w
i
θ
2
+ ν

(28)
em que λ
−i
denota todos os elementos de λ, exceto λ
i
. Verifica-se a partir de
(28) que cada λ
i
tem distribuição condicional proporcional a uma gamma com
parâmetros dados abaixo.
λ
i
|λ
−i
, L
∗
, θ
∗
, ν ∼ G
ν + 1
2
,
1
2
L
∗
i
− w
i
θ
2
+ ν

. (29)
A distribuição condicional para L
i
resultará em uma normal conforme a
expressão:
L
∗
i
|λ
i
, θ
∗
, ν ∼ N
w
i
θ
∗
,
1
λ
i
, γ
∗
k−1
< L
∗
i
< γ
∗
k
, k = 1, 2, 3, ..., K. (30)
A distribuição condicional completa para ν não tem forma fechada e assim,
torna-se necessário o uso do algoritmo Metropolis-Hastings para a amostragem
19
desta distribuição
p(ν|θ
∗
, L
∗
, γ
∗
, λ) ∝
ν
2
(ν/2)
Γ (ν/2)
n
n
i=1
λ
ν
2
−1
i
exp
−
ν
2
λ
i
1
(1 + ν)
2
(31)
As distribuições condicionais para os parâmetros σ
2
u
e γ
k
são, respectiva-
mente, dadas pelas expressões (15) e (25).
Neste caso, a amostragem Gibbs é realizada por meio das condicionais
completas dadas em (27), (15), (18), (29), (31) e (25), não necessariamente nesta
ordem.
Dada a dificuldade de se obter um valor adequado para a variância da dis-
tribuição geradora de candidatos para (γ
∗
k
), na metodologia proposta por Cowles
(1996) e com o intuito de melhorar o processo de convergência, Nandran & Chen
(1996) apresentam a reparametrização dada por:
δ = 1/γ
∗
K−1
, γ
∗∗
k
= δγ
∗
k
, k = 0, 1, 2, ..., K, θ
∗∗
= δθ
∗
e L
∗∗
= δL
∗
. (32)
A amostragem de γ
∗∗
e L
∗∗
é como em Cowles (1996). Com esta repara-
metrição, os parâmetros de limiar a serem estimados (γ
∗
k
) ficam limitados entre 0
e 1. Além disso, utilizando o teorema do valor médio, verifica-se que a distribui-
ção de (γ
∗
|θ
∗
, Y ) é dada de forma proporcional ao produto de uma distribuição
normal e uma ditribuição Dirichlet. A limitação dos parâmetros de limiar entre
0 e 1, motiva à utilização desta última como geradora de candidatos para estes
parâmetros, ainda que indiretamente. Outro atrativo desta forma de gerar candida-
tos é que com esta distribuição o vetor dos valores propostos é gerado de uma só
vez e não depende dos demais parâmetros. Nesta abordagem, a função de ligação
considerada para (2) é a distribuição normal acumulada e um modelo fixo é ado-
tado para L
∗∗
. Um dos nossos propósitos é adaptar a metodologia de Nandran &
20
TABELA 2: Formas de análise encontradas na literatura para o problema de dados
categorizados em escalas subjetivas.
Função geradora Distribuição Forma de análise
de candidatos da variável Modelos fixos Modelos com efeitos
para γ latente aleatórios
Condicionais Normal Sörensem et al. (1995) Sörensem et al. (1995)
completas t de student Albert & Chib (1993)
Normal Normal Cowles (1996) Kizilkaya et al.(2003)
t de Student Kizilkaya et al.(2003)
Dirichlet Normal Nandran & Chen (1996)
t de Student
Chen (1996) para modelos mistos, com distribuições normal e t de Student para a
variável latente.
Na Tabela 2 estão indicados os algoritmos, as distribuições adotadas para
a variável latente, o tipo de modelo (fixo ou aleatório), bem como, a autoria en-
contrados na literatura. As células em branco na coluna, que se refere à forma de
análise por meio de modelos com efeitos aleatórios, compreendem às metodolo-
gias que serão desenvolvidas neste trabalho. Nesta Tabela está indicado também
que a amostragem Gibbs é utilizada nos algoritmos desenvolvidos por Sörensen
et al. (1995) e Albert & Chib (1993), e nos outros dois a correspondente distri-
buição geradora de candidatos, pois, a distribuição conjunta da variável latente e
dos parâmetros threshold é escrita de forma fatorada e faz-se necessário o uso do
algoritmo Metropolis-Hastings.
As formas de análises apresentadas na Tabela 2 podem ser classificadas
em três grupos quanto à amostragem. O primeiro utiliza a amostragem Gibbs, o
segundo a amostragem por meio do algoritmo Metropolis-Hastings (MH) e distri-
buição normal como geradora de candidatos e o último, a amostragem utilizando
o algoritmo MH e distribuição Dirichlet como geradora de candidatos. Estas três
formas de amostragem foram originalmente desenvolvidas, respectivamente, por
21
Albert & Chib (1993), Cowles (1996) e Nandram & Chen (1996), como já men-
cionado. As outras autorias indicadas na Tabela 2 são adaptações quanto à distri-
buição adotada para a variável latente e modelos com efeitos aleatórios.
2.4 Comparação de Modelos
A comparação entre modelos pode ser realizada empregando o fator de
Bayes (Gelman et al., 2003), o qual de forma geral é dado por:
B
ij
=
p(y|M
i
)
p(y|M
j
)
=
p(y|θ
i
, M
i
)p
i
(θ
i
|M
i
)dθ
i
p(y|θ
j
, M
j
)p
j
(θ
j
|M
j
)dθ
j
(33)
em que θ
i
e θ
j
são os vetores de parâmetros dos modelos M
i
e M
j
, respectiva-
mente, com M
i
e M
j
dois modelos concorrentes. Se o fator (B
ij
) for maior do
que 1 há evidências de que o modelo i seja melhor do que o modelo j.
Porém, nem sempre é possível obter as integrais indicadas em (33), alter-
nativamente o método da amostragem por importância pode ser considerado para
obter uma aproximação da marginal dos dados e, neste caso, a comparação dos
modelos se torna mais fácil, pois basta obter uma estimativa para a média ou para
a média harmônica da verossimilhança caso a função de importância (g(θ
i
)) seja,
respectivamente, a distribuição a priori e a distribuição a posteriori dos parâme-
tros.
As duas formas de estimar a marginal dos dados, apresentadas a seguir, são
descritas segundo Sörensen & Gianola (2004). A distribuição marginal dos dados
22
pode ser expressa como
p(y|M
i
) =
p(y|θ
i
, M
i
)p(θ
i
|M
i
)dθ
i
p(θ
i
|M
i
)dθ
i
=
p (y|θ
i
, M
i
)
p(θ
i
|M
i
)
g(θ
i
)
g (θ
i
) dθ
i
p(θ
i
|M
i
)
g(θ
i
)
g (θ
i
) dθ
i
(34)
o denominador em (34) reescrito, como a seguir, pode ser interpretado como a
esperança de p(θ
i
|M
i
)/g(θ
i
).
p(θ
i
|M
i
)
g(θ
i
)
g(θ
i
)dθ
i
= lim
m→∞
1
m
m
j=1
p
θ
[j]
i
|M
i
g
θ
[j]
i
(35)
em que p
θ
[j]
i
|M
i
é a densidade a priori para o modelo i avaliada no j-ésimo
valor amostrado para θ.
O numerador em (34) de forma análoga ao denominador pode ser escrito
como:
p(y|θ
i
, M
i
)
p(θ
i
|M
i
)
g(θ
i
)
g(θ
i
)dθ
i
= lim
m→∞
1
m
m
j=1
p
y|θ
[j]
i
, M
i
p
θ
[j]
i
|M
i
g
θ
[j]
i
(36)
em que p
y|θ
[j]
i
, M
i
é avaliada no j-ésimo valor obtido da distribuição de impor-
tância.
Assim, fazendo ω = p
θ
[j]
i
|M
i
/g
θ
[j]
i
, para m grande a distribuição
marginal dos dados é estimada por meio da razão entre (36) e (35) como a seguir:
p(y|M
i
) =
m
j=1
ω
[j]
i
p
y|θ
[j]
i
, M
i
m
i
ω
[j]
i
(37)
23
TABELA 3: Escala para interpretação do fator de Bayes.
B
ij
Evidência favorável a M
i
< 1 negativa (favorável a M
j
)
de 1 a 3 duvidosa
de 3 a 10 substancial
de 10 a 30 forte
de 30 a 100 muito forte
> 100 decisiva
Fonte: Jeffreys (1961).
Se a função de importância for a distribuição a priori, então cada um dos
pesos (ω
i
) será igual a 1 e a distribuição marginal dos dados (37) resultará em:
p(y|M
i
) =
1
m
m
j=1
p
y|θ
[j]
i
, M
i
, (38)
com θ
[j]
i
retirado da distribuição a priori.
A vantagem do uso da distribuição a priori como função de importância é
a facilidade de obtenção de uma amostra de θ e a desvantagem é que, em geral,
com esses parâmetros, encontra-se valores muito pequenos para a p(y|θ).
Por outro lado, se a função de importância for a distribuição a posteriori,
obtém-se maior probabilidade para os dados e, além disso, não há necessidade
do conhecimento da forma desta distribuição, conforme a expressão em (40). A
desvantagem é a instabilidade numérica para o estimador da distibuição marginal
dos dados (p(y|M
i
)), mas esse problema pode ser contornado utilizando a escala
logarítmica.
ω =
p (θ
i
|M
i
)
p (θ
i
|y, M
i
)
=
p (θ
i
|M
i
)
p(y|θ
i
,M
i
)p(θ
i
|M
i
)
p(y|M
i
)
=
p (y|M
i
)
p (y|θ
i
, M
i
)
(39)
24
e substituindo a expressão à direita da última igualdade em (37) obtêm-se:
p(y|M
i
) =
m
j=1
p(y|M
i
)
p
y|θ
[j]
i
,M
i
p
y|θ
[j]
i
, M
i
m
j=1
p(y|M
i
)
p
y|θ
[j]
i
,M
i
=
m
m
j=1
1
p
y|θ
[j]
i
,M
i
=
1
m
m
j=1
1
p
y|θ
[j]
i
, M
i
−1
. (40)
O estimador da distribuição marginal dos dados em (40) é a média harmô-
nica dos valores da verossimilhança calculada em cada ciclo j do processo MCMC.
A interpretação do fator de Bayes segundo a escala dada por Jeffreys (1961)
é apresentada na Tabela 3. Com esta escala é possível visualizar, dentre dois mo-
delos em comparação, a intensidade com a qual um é superior ao outro.
25
3 MATERIAL E MÉTODOS
3.1 Material
Os dados (Tabela 4) utilizados para a ilustração da metodologia proposta
neste trabalho foram obtidos pelo Centro Federal de Educação Tecnológica de Rio
Verde (CEFET-RV) durante a realização da Exposição Agropecuária deste muni-
cípio no ano de 2006. A variável resposta consistiu de notas atribuídas por 36
provadores a amostras de banana da terra desidratada. Os tratamentos avaliados
nesta análise foram três diferentes concentrações de sacarose utilizadas no pro-
cesso de desidratação osmótico, sendo T1 - solução a 30%, T2 - solução a 40% e
T3 - solução a 50% de sacarose.
Os provadores incluíam crianças e adultos de ambos os sexos, não treina-
dos. Empregou-se a escala hedônica de nove pontos, cujo ponto 1 correspondeu
a “desgostei muitíssimo” e o 9 a “gostei muitíssimo” (Dutcosky, 1996). A ca-
racterística analisada foi a cor da banana da terra desidratada. O delineamento
experimental adotado foi em blocos ao acaso em que cada provador avaliou os três
tratamentos e, assim, este constituía um bloco completo.
3.2 Métodos
Dentre as distribuições que podem ser empregadas para explicar o fenô-
meno estudado, neste trabalho focou-se na comparação entre a t de Student e a
normal.
Seguindo recomendação de Kizilkaya et al. (2002) as categorias foram
reagrupadas de tal forma que em cada uma tivesse no mínimo cinco observações.
O esquema de reorganização está descrito na Tabela 5.
26
TABELA 4: Notas atribuídas por 36 provadores à cor da banana da terra desidra-
tada utilizando três concentrações de sacarose.
Provador % de sacarose Provador % de sacarose
30 40 50 30 40 50
1 7 1 2 19 8 7 8
2 8 8 8 20 9 9 9
3 9 6 5 21 7 6 7
4 7 3 8 22 7 8 8
5 7 8 8 23 9 9 9
6 9 9 9 24 8 4 9
7 9 8 8 25 9 8 8
8 9 9 9 26 9 8 8
9 9 9 6 27 9 9 9
10 8 7 7 28 8 8 8
11 8 9 8 29 8 7 6
12 8 5 9 30 9 9 9
13 9 8 9 31 8 8 8
14 8 6 5 32 8 8 8
15 9 9 9 33 8 8 9
16 9 6 8 34 8 9 6
17 8 9 7 35 5 5 5
18 8 8 8 36 7 9 8
Fonte: Dados fornecidos pelo Centro Federal de Educação Tecnológica de Rio
Verde - CEFET-RV.
TABELA 5: Esquema de reclassificação dos dados em cinco categorias.
Categorias E. H. de 9 Pontos* E. H. de 5 Pontos
Desgostei muitíssimo 1
Desgostei muito 2
Desgostei 3 1
Desgostei moderadamente 4
Nem gostei nem desgostei 5
Gostei moderadamente 6 2
Gostei 7 3
Gostei muito 8 4
Gostei muitíssimo 9 5
(*)Fonte: Dutcosky (1996).
27
A análise de variância (ANAVA)usual foi realizada para os dados sem
transformação, com transformação logarítmica e com transformação de Box-Cox
(Box & Cox, 1964). Também, foi apresentada a análise de variância para o modelo
de regressão linear. A pressuposição de normalidade dos resíduos do modelo em
cada transformação, foi estudada por meio do gráfico dos resíduos em função dos
quantis teóricos (Normal Q-Q Plot) e, por meio do teste de Shapiro-Wilk (Shapiro
& Wilk, 1965).
Na análise da variância, a variável resposta Y , foi modelada como segue:
Y
ij
∼ N
β
i
+ u
j
, σ
2
,
em que Y
ij
é a observação ou nota atribuída ao tratamento i (i=1, 2, 3), pelo pro-
vador (bloco) j (j = 1, 2, ··· , 36); β
i
é a média do tratamento i, u
j
o efeito do
bloco j e σ
2
a variância residual.
Os elementos de Y são observados em uma das K categorias ordenadas e
fazendo p
ik
= P [Y
i
= k] pode-se definir as probabilidades acumuladas por meio
da expressão η
ik
=
K−1
k=1
p
ik
. Segundo McCullagh (1980) um modelo para p
ik
é dado por η
ik
= Φ(γ
k
− x
i
θ), com i = 1, 2, ··· , n e k = 1, , 2, ··· , K − 1
e generalizando, η
ik
= G(γ
k
− x
i
θ). Portanto, verifica-se que a análise pode
ser realizada por meio dos modelos lineares generalizados (GLM) cuja função de
ligação é uma distribuição de probabilidade acumulada. Neste trabalho, foram
atribuídas as distribuições normal e t de Student acumuladas para G.
Para a implementação do processo de Monte Carlo via Cadeias de Markov
(MCMC), no caso em que se assumiu distribuição normal para a variável latente
28
L, o modelo para Y foi escrito como:
p(Y
i
= k|θ, γ) = Φ
γ
k
− w
i
θ
δ
2
− Φ
γ
k−1
− w
i
θ
δ
2
em que i = 1, 2, ··· , n, com n representado o número de observações e k =
1, 2, ··· , K é o índice para categorias, Φ(·) a distribuição normal padrão acumu-
lada e θ
= c(β
, u
). Introduzindo a variável latente L o modelo foi reescrito
como segue:
p(Y
i
= k|L
i
, γ) =
K
k=1
I
[γ
k−1
;γ
k
]
(L
i
)I(Y
i
= k)
com
L
i
|θ, δ
2
∼ N(w
i
θ, δ
2
)
Quando assumiu-se disribuição t de Student para L, a especificação do
modelo foi como a seguir:
p(Y
i
= k|θ, γ) = F
ν
γ
k
− w
i
θ
δ
2
− F
ν
γ
k−1
− w
i
θ
δ
2
em que F
ν
representa a distribuição t de Student acumulada. O modelo em dois
estágios resultou em:
p(Y
i
= k|L
i
, γ) =
K
k=1
I
[γ
k−1
;γ
k
]
(L
i
)I(Y
i
= k)
com
29
L
i
|θ, δ
2
∼ N(w
i
θ,
δ
2
λ
i
)
e
p(λ
i
|ν) =
ν
2
ν
2
Γ
ν
2
λ
ν
2
−1
i
exp
−
λ
i
2
ν
(41)
Nos algoritmos AC e MC foi assumido, sem perda de generalidade, que
δ
2
= 1.
Foi considerada informação a priori vaga para β por meio da distribuição,
p(β) ∝ exp
−
0, 001
2
β
β
A especificação da priori para os efeitos aleatórios u foi feita por meio de
uma distribuição normal multivariada como segue:
u|A ∼ N (∅, Aσ
2
u
),
em que ∅é um vetor nulo, de dimensão q×1. Em todas as abordagens consideradas
neste trabalho, a matriz A é uma matriz identidade q × q por assumir que em
experimento com blocos aleatorizados estes são não correlacionados.
As prioris para ν, δ
2
e σ
2
u
foram respectivamente:
p(ν) ∝
1
(1 + ν)
2
,
p
σ
2
u
∝
σ
2
u
−(20+1)
exp
−5
σ
2
u
, (42)
e
30
p
σ
2
u
∝
σ
2
u
−(3+1)
exp
−5
σ
2
u
(43)
Considerou-se as mesmas prioris em todos os algoritmos meramente para
fins de padronização eliminando, assim, a possibilidade de um algoritmo ser pri-
vilegiado por prioris que contribuíssem para a convergência mais rápida ou preju-
dicado por prioris que produzissem efeito contrário.
Os Três tratamentos do exemplo são níveis quantitativos igualmente es-
paçados do fator concentração de sacarose. Assim, o modelo de regressão linear
simples pode ser ajustado com o modelo quadrático permitindo testar a falta de
ajustedo modelo linear.
As médias de tratamentos podem ser estruturadas em dois contrastes orto-
gonais referentes aos efeitos de regressão linear e quadrático, quais sejam:
• C
1
β =
β
3
− β
1
2
, associado à regressão linear para efeito de porcentagem de
sacarose; e
• C
2
β = β
2
−
β
3
+ β
1
2
, associado à regressão quadrática.
Tais contrastes podem ser diretamente calculados na amostra da distribui-
ção a posteriori conjunta.
O coeficiente de correlação intraclasse (ρ),
ρ =
σ
2
u
σ
2
u
+ δ
2
foi calculado para avaliar a semelhança entre as notas dos diferentes avaliadores.
As metodologias propostas por Sörensen et al. (1995) (seção 2.1) e por Ki-
zilkaya et al. (2003) adaptada de Cowles (1996) (seções 2.2 e 2.3) foram aplicadas
aos dados descritos anteriormente. Foram desenvolvidas, também, adaptações para
modelos mistos, a partir da metodologia de Albert & Chib (1993) (seção 4.1.1) e
31
para modelos mistos com funções de ligação normal (seção 4.1.2) e t de Student
(seção 4.1.2) acumuladas a partir do algoritmo de Nandram & Chen (1996).
Os algoritmos foram implementados em R (R Development Core Team,
2007) e permitem analisar quaisquer estruturas de covariância para os efeitos ale-
atórios.
Por meio da distribuição conjunta a posteriori em cada uma da situações
obteve-se cadeias com 10000 iterações para cada parâmetro. Como um dos objeti-
vos deste trabalho foi a verificação da convergência nos três algorítmos (AC, MC
e NC) e nas duas funções de ligação utilizadas, para este fim, não foi considerada
uma amostra válida para a realização de inferências e sim, uma cadeia de tamanho
suficiente para o estudo da convergência em cada caso. Para a aplicação do critério
de Gelman & Rubin foram consideradas duas cadeias para cada parâmetro.
Foram comparadas também a modelagem por meio da distribuição t de
Student e por meio da distribuição normal. Nestes casos, foi necessário o uso
da amostra válida, pois, além da velocidade de convergência dos algoritmos, o
interesse também consistiu na verificação da melhor modelagem (distribuição t de
Student ou normal). A escolha do melhor modelo foi realizada em cada um dos
algoritmos, por meio do fator de Bayes, estimado por meio da média harmônica
comforme a expressão em (40).
Para a realização das inferências foi considerada uma amostra final de ta-
manho 5000, obtida a partir de uma cadeia com 255000 iterações, em que as 5000
primeiras foram descartadas para a eliminação do efeito do valor inicial arbitrário
e armazenado um valor a cada 50 para que estes fossem não correlacionados.
Essas cadeias foram utilizadas também para a estimação dos contrastes
anteriormente citados.
A verificação da convergência foi realizada por meio de gráficos para as
32
cadeias dos parâmetros, pelos critérios de Raftery & Lewis (1992) e de Gelman &
Rubin (1992).
Foi calculada a autocorrelação para todos os parâmetros no processo amos-
tral. O cálculo dos valores de autocorrelações para diferentes intervalos (autocor-
relação de lag k) entre os dados das cadeias de Markov foi feito como se segue:
r
k
=
E[(X
t
− X)(X
t+k
− X)]
E[(X
t
− X)
2
]E[(X
t+k
− X)
2
]
em que X
t
é o valor amostrado na iteração t, X
t+k
o valor amostrado na iteração
t + k e X é a média de todos os valores amostrados.
Estão apresentadas autocorrelações calculadas para os parâmetros corres-
pondentes à média do tratamento e ao parâmetro threshold que apresentaram maior
problema de autocorrelação em cada algoritmo.
33
4 RESULTADOS E DISCUSSÃO
4.1 Resultados metodológicos
Nesta seção são apresentados os desenvolvimentos das três adaptações
aplicadas ao modelo misto, realizadas a partir das metodologias de Nandran &
Chen (1996), com função de ligação normal e t de Student acumuladas e de Albert
& Chib (1993) com função de ligação t de Student.
Na seção 4.2 será apresentado um exemplo de aplicação dos modelos de
Albert & Chib (1993), Cowles (1996) e Nandran & Chen (1996), bem como suas
adaptações na análise sensorial de alimentos. O estudo de convergência é realizado
em cada caso.
4.1.1 Algoritmo adaptado de AC (AC-t)
A metodologia de Albert & Chib (1993) é adaptada neste caso, com a
introdução de parâmetros de efeitos aleatórios na modelagem de L
∗
.
As distribuições condicionais completas para cada parâmetro foram obti-
das a partir da equação (24) e são as que seguem:
θ
∗
|L
∗
, λ, ν ∼ N
M
−1
W
R
−1
L
∗
, M
−1
(44)
em que M = W
R
−1
W + V
−1
e R é uma matriz diagonal com seus elementos
dados por λ
−1
i
i = 1, 2, ··· , n e,
p(λ
i
|λ
−i
, L
∗
, θ, ν) ∝ λ
(
ν+1
2
)
−1
i
exp
−
λ
i
2
L
∗
i
− w
i
θ
2
+ ν

(45)
a qual é o núcleo de uma distribuição Gama com parâmetros dados a seguir e λ
−i
,
34
denota todos os elementos de λ exceto λ
i
.
λ
i
|λ
−i
, L
∗
, ν, θ ∼ G
ν + 1
2
,
1
2
L
∗
i
− w
i
θ
2
+ ν

(46)
e
L
∗
i
|L
∗
−i
, θ, ν, λ
i
∼ N
w
i
θ
∗
,
1
λ
i
, γ
∗
k−1
< l
∗
i
< γ
∗
k
, k = 1, 2, 3, ..., K. (47)
A distribuição condicional completa a posteriori para ν não tem forma
fechada e, assim, torna-se necessário o uso do algoritmo Metropolis-Hastings para
a amostragem desta distribuição:
p(ν|θ
∗
, L
∗
, γ
∗
) ∝
ν
2
(ν/2)
Γ (ν/2)
n
n
i=1
λ
ν
2
−1
i
exp
−
ν
2
λ
i
(48)
As distribuições condicionais para os parâmetros σ
2
u
e γ
k
são, respectiva-
mente, dadas pelas expressões (15) e (11).
A implementação do Amostrador de Gibbs pode ser feita amostrando se-
guidamente das distribuições condicionais completas, dadas em (44), (15), (46),
(48), (47) e (11).
4.1.2 Algoritmo adaptado de NC - (NC-normal)
Nesta seção é apresentada a proposta de Nandram & Chen (1996) adaptada
para o modelo misto, isto é, a variável latente é modelada em função de efeitos
fixos e aleatórios. Estes autores propuseram a reparametrização conforme (49) e
utilizaram a distribuição Dirichlet como candidata para os parâmetros de limiar,
δ = 1/γ
∗
K−1
, γ
∗∗
k
= δγ
∗
k
, k = 0, 1, 2, ..., K, θ
∗∗
= δθ
∗
e L
∗∗
= δL
∗
. (49)
35
Desta forma, os parâmetros threshold continuam dividindo a reta em K in-
tervalos e destes parâmetros, agora, existem apenas K −3 desconhecidos. Observe
que com essa reparametrização, para problemas com três categorias, não existirá
nenhum parâmetro de limiar a ser estimado.
O Jacobiano da transformação é dado por
δ
2
−
1
2
(n+m+K)
, em que n é o
número de observações, m o número de variáveis explicativas e K a quantidade de
categorias. Considerando uma distribuição a priori para δ
2
, a distribuição conjunta
a posteriori em (8) fica sendo:
p
θ
∗∗
, γ
∗∗
k
, δ
2
, L
∗∗
|y
∝
n
i=1
φ
L
∗∗
i
, w
i
θ
∗∗
, δ
2
I
[γ
∗∗
k
;γ
∗∗
k+1
)
(L
∗∗
i
)
× φ
θ
∗∗
; 0, δ
2
V
δ
2
−
K
2
p
σ
2
u
p
δ
2
(50)
A distribuição condicional completa a posteriori para θ
∗∗
é,
θ
∗∗
|L
∗∗
, δ
2
, Y, σ
2
u
∼ N
B
−1
W
L
∗∗
, δ
2
B
−1
. (51)
com B, neste caso, dado por: B = δ
2
V
−1
+ W
W
Se p
δ
2
na expressão anterior for uma Gamma Inversa,
p
δ
2
∝
δ
2
−(c+1)
exp
−
d
δ
2
(52)
em que, c e d são os hiperparâmetros desta priori, então, a distribuição condicional
completa a posteriori para δ
2
será:
δ
2
|θ
∗∗
, L
∗∗
, Y, σ
2
u
∼ GI (a
δ
, b
δ
) (53)
36
com
a
δ
=
(n + m + K + 2c)
2
e
b
δ
=
(L
∗∗
− W θ
∗∗
)
(L
∗∗
− W θ
∗∗
) + θ
∗∗
V
−1
θ
∗∗
+ 2d
2
A variância de efeitos aleatórios (σ
2
u
) tem distribuição condicional como
em (15).
Como já mencionado, a distribuição condicional de (γ
∗∗
, L
∗∗
|θ
∗∗
, δ
2
, Y )
pode ser escrita como o produto das distribuições condicionais (γ
∗∗
|θ
∗∗
, δ
2
, Y ) e
(L
∗∗
|γ
∗∗
, θ
∗∗
, δ
2
, Y ). E com a reparametrização em (49) a distribuição condicio-
nal para (γ
∗∗
|θ
∗∗
, δ
2
, Y ) é
π
γ
∗∗
|θ
∗∗
, δ
2
, Y
∝
Y
i
=2
Φ
γ
∗∗
2
− w
i
θ
∗∗
δ
− Φ
−w
i
θ
∗∗
δ

×
Y
i
=3
Φ
γ
∗∗
3
− w
i
θ
∗∗
δ
− Φ
γ
∗∗
2
− w
i
θ
∗∗
δ

···
··· ×
Y
i
=K−1
Φ
1 − w
i
θ
∗∗
δ
− Φ
γ
∗∗
K−2
− w
i
θ
∗∗
δ

(54)
e para
L
∗∗
|γ
∗∗
, θ
∗∗
, δ
2
, Y
L
∗∗
|γ
∗∗
, θ
∗∗
, δ
2
, Y
i
= k ∼ N
W θ
∗∗
, δ
2
. (55)
Observe que para a amostragem de γ
∗∗
é necessário o uso do algoritmo
37
Metropolis-Hastings, pois este parâmetro não tem uma distribuição com forma na
qual se possa aplicar o método de Monte Carlo. Cowles (1996) propôs a distribui-
ção normal, para a geração de candidatos para um determinado threshold, truncada
por outros dois adjacentes e Nandram & Chen (1996) propuseram a distribuição
Dirichlet, a qual é uma distribuição conjugada da multinomial.
Uma motivação para o uso da distribuição Dirichlet é que os threshold, na
forma reparametrizada, os valores de (γ
∗∗
k
) ficam no intervalo de zero a um, isto é,
a partição em (2) resultará em:
(γ
∗∗
0
, γ
∗∗
1
); [γ
∗∗
1
, γ
∗∗
2
); ··· ; [γ
∗∗
K−1
, γ
∗∗
K
);
com γ
∗∗
1
= 0 e γ
∗∗
K−1
= 1 e assim, pode se construir uma variável p definida da
seguinte forma:
p
k−1
= γ
∗∗
k
− γ
∗∗
k−1
, k = 2, ..., K − 1;
em que
p = (p
1
, p
2
, ..., p
K−2
)
, p
k
≥ 0, k = 1, 2, ..., K − 2 e
K−2
k=1
p
k
= 1
De acordo com o teorema do valor médio,
Φ
γ
∗∗
k
− w
i
θ
∗∗
δ
− Φ
γ
∗∗
k−1
− w
i
θ
∗∗
δ
=
1
δ
φ
ξ
k−1
− w
i
θ
∗∗
δ
p
k−1
(56)
em que ξ
k−1
∈
γ
∗∗
k
; γ
∗∗
k−1
, k = 2, 3, ..., K − 1, e φ (·) é a função densidade
normal padrão.
Então, a partir de (56) tem-se:
π
γ
∗∗
|θ
∗∗
, δ
2
, Y
∝ h
1
(ξ) h
2
(p) (57)
38
e na expressão anterior, h
1
e h
2
são, respectivamente:
h
1
(ξ) =
K−2
k=1
n
k
i=1
φ
ξ
k
− w
i
θ
∗∗
δ
e
h
2
(p) =
K−2
k=1
p
n
k+1
k
(58)
Nesta adaptação, a distribuição Dirichlet também é usada como geradora
de candidatos pois, a reparametrização dos parâmetros threshold é como desen-
volvido por Nandram & Chen (1996).
Pode-se observar que a distribuição condicional de γ
∗∗
é dada pelo produto
de h
1
por h
2
em que h
2
é o núcleo de uma distribuição distribuição Dirichlet com
parâmetros n = (n
2
+1, n
3
+ 1, ...,n
k−1
+ 1)
e h
2
não depende de θ
∗∗
e nem de
δ.
Com os novos valores propostos para p gerados a partir de h
2
a construção
dos novos valores de γ
∗∗
k
é como segue:
γ
∗∗
k,j
=
k−1
i=1
p
i,j
, k = 2, 3, ..., K − 2.
A probabilidade de aceitação do novo vetor
γ
∗∗
j
=
γ
∗∗
2,j
, γ
∗∗
3,j
, ..., γ
∗∗
K−2,j
é o mínimo entre os valores um e α em que
α = ω
γ
∗∗
j
, p
j
/ω
γ
∗∗(j−1)
, p
j−1
, (59)
e j representa a j-ésima iteração, do algoritmo. Em (59), a forma geral da dis-
tribuição de ω (γ
∗∗
, p) é dada por ω (γ
∗∗
, p) = p
γ
∗∗
|θ
∗∗
, δ
2
/p
p|n, θ
∗∗
, δ
2
.
39
A expressão de p(γ
∗∗
|θ
∗∗
, δ
2
) é dada em (54) e p(p|n, θ
∗∗
, δ
2
) é a distribuição
Dirichlet,
p(p) =
1
Z(n)
K−2
k=1
p
n
k+1
−1
k
(60)
em que p
1
, ··· , p
K−2
≥ 0;
K−2
k=1
p
k
= 1 e n
2
, ··· , n
K−1
> 0. Na expressão
acima, Z(n) é a constante normalizadora dada por
Z(n) =
K−2
k=1
Γ(n
k+1
)
Γ(
K−2
k=1
n
k+1
)
.
Como o interesse é na distribuição conjunta de L
∗∗
e γ
∗∗
, a amostragem é feita
para p e, a partir dos elementos deste vetor constrói-se γ
∗∗
k
, e caso o novo vetor
de parâmetros threshold seja aceito atualiza-se os valores de L
∗∗
. Caso contrário
continua com a amostra anterior da variável latente.
Para a implementação deste algoritmo, a amostragem é realizada a partir
das expressões em (51), (15), (53), (60) e (55).
4.1.3 Algoritmo adaptado de NC - (NC-t)
Na forma reparametrizada, dada em (49), ou seja, de acordo com a metodo-
logia apresentada por Nandram & Chen (1996) (algoritmo NC) e com a utilização
da distribuição t de Student acumulada como função de ligação, a verossimilhança
fica sendo:
P (Y
i
= k|θ
∗∗
, δ, ν, γ
∗∗
) = F
ν
γ
∗∗
k
− w
i
θ
δ
− F
ν
γ
∗∗
k−1
− w
i
θ
δ
, (61)
k = 1, 2, ..., K
em que F
ν
é a distribuição t-Student acumulada com ν graus de liberdade.
40
Como já visto, a distribuição t de Student pode ser escrita em dois estágios,
os quais são dados agora, como segue:
L
∗∗
i
|θ
∗∗
, δ
2
, λ
i
∼ N
w
i
θ
∗∗
,
δ
2
λ
i
, i = 1, 2, ..., n (62)
λ
i
|ν ∼ Gamma
ν
2
,
ν
2
. (63)
O modelo amostral em (61) será dado por:
P (Y
i
= k|θ
∗∗
, δ, λ
i
, ν, γ
∗∗
) = Φ
γ
∗∗
k
− w
i
θ
δ
√
λ
i
− Φ
γ
∗∗
k−1
− w
i
θ
δ
√
λ
i
, (64)
k = 1, 2, ..., K
Para a obtenção da distribuição condicional completa própria para ν, Kizilkaya et
al. (2003) propuseram a priori p(ν) = 1/(1 + ν)
2
e admitindo priori vaga para
θ
∗∗
a posteriori conjunta de todos os parâmetros e a variável auxiliar ou latente é:
p(θ
∗∗
, γ
∗∗
k
, δ
2
, L
∗∗
, λ|y) ∝
∝ (δ
2
)
−
K
2
n
i=1
φ
L
∗∗
i
, w
i
θ
∗∗
,
δ
2
λ
i
I
[γ
∗∗
k
,γ
∗∗
k+1
)
(L
∗∗
i
)
×
n
i=1
λ
(
ν
2
)
−1
i
exp
−
λ
i
2
ν
φ
θ
∗∗
, 0, δ
2
V
1
(1 + ν)
2
× p
σ
2
u
p
δ
2
p(ν) (65)
em que λ = {λ
i
}
n
i=1
, V é como na seção 2.1, e as prioris para σ
2
u
e δ
2
são dadas
como em (7) e (52), respectivamente.
As distribuiões condicionais completas a posteriori, para cada parâmetro
são dadas a seguir.
41
A distibuição condicional para θ, resulta em:
θ
∗∗
|L
∗∗
, δ
2
, ν, λ ∼ N
M
−1
W
R
−1
L
∗∗
, δ
2
M
−1
(66)
em que M = W
R
−1
W + δ
2
V
−1
.
O parâmetro λ
i
tem distribuição condicional conforme expressão a seguir:
p(λ
i
|λ
−i
, L
∗∗
, θ, ν) ∝ λ
(
ν+1
2
)
−1
i
exp
−
λ
i
2
L
∗∗
i
− w
i
θ
2
+ ν

(67)
com λ
−i
, denotando todos os elementos de λ exceto λ
i
. A distribuição acima é
proporcional a uma Gama com parâmetros (ν + 1)/2 e
L
∗∗
i
− w
i
θ
2
+ ν, isto é,
λ
i
|λ
−i
, L
∗∗
i
, θ, ν, δ
2
∼ G
ν + 1
2
,
1
2
L
∗∗
i
− w
i
θ
2
+ ν

(68)
e
L
∗∗
i
|λi, θ, ν, δ
2
∼ N
w
i
θ
∗∗
,
δ
2
λ
i
, γ
k−1
< l
∗∗
i
< γ
k
, k = 1, 2, 3, ..., K (69)
A distribuição condicional completa para ν não tem forma fechada e, assim, torna-
se necessário o uso do algoritmo Metropolis-Hastings para a amostragem desta
distribuição.
p(ν|θ
∗∗
, L
∗∗
, γ, δ
2
) ∝
ν
2
(ν/2)
Γ (ν/2)
n
n
i=1
λ
ν
2
−1
i
exp
−
ν
2
λ
i
1
(1 + ν)
2
. (70)
A distribuição condicional completa para δ
2
tem a mesma forma como em (53).
Como para o caso do modelo normal, a distribuição conjunta de L
∗∗
e
γ
∗∗
pode ser escrita por meio do produto das condicionais (γ
∗∗
|θ
∗∗
, δ
2
, λ, Y ) e
(L
∗∗
|γ
∗∗
, θ
∗∗
, δ
2
, λ, Y ).
42
A distribuição condicional completa para a variável auxiliar é como em
(69). Para γ
∗∗
a distribuição condicional completa é dada por:
p(γ
∗∗
|θ
∗∗
, δ
2
, λ, Y ) ∝
Y
i
=2
Φ
γ
∗∗
2
− w
i
θ
∗∗
δ
√
λ
i
− Φ
−w
i
θ
∗∗
δ
√
λ
i

×
Y
i
=3
Φ
γ
∗∗
3
− w
i
θ
∗∗
δ
√
λ
i
− Φ
γ
∗∗
2
− w
i
θ
∗∗
δ
√
λ
i

···
Y
i
=K−1
Φ
1 − w
i
θ
∗∗
δ
√
λ
i
− Φ
γ
∗∗
K−2
− w
i
θ
∗∗
δ
√
λ
i

(71)
E neste caso também, de acordo com o teorema do valor médio, têm-se:
Φ
γ
∗∗
k
− x
i
θ
∗∗
δ
√
λ
i
− Φ
γ
∗∗
k−1
− x
i
θ
∗∗
δ
√
λ
i
=
1
δ
√
λ
i
φ
ξ
k−1
− x
i
θ
∗∗
δ
√
λ
i
p
k−1
(72)
em que ξ
k−1
∈ (γ
∗∗
k
; γ
∗∗
k−1
), k = 2, 3, ..., K − 1, e φ(·) é a função densidade
normal padrão.
E a partir de (72)
π(γ
∗∗
|θ
∗∗
, δ
2
, Y ) ∝ h
3
(ξ)h
4
(p) (73)
e em (71), h
3
e h
4
são, respectivamente:
h
3
(ξ) =
K−2
k=1
n
k
i=1
ξ
k
− w
i
θ
∗∗
δ
√
λ
i
e
h
4
(p) =
K−2
k=1
p
n
k+1
k
(74)
No caso de adotar a distribuição t de Student para a variável latente, os
43
procedimentos para gerar os candidatos, a construção dos threshold e a probabili-
dade de aceitação (59) são idênticos ao caso do modelo normal, com exceção do
parâmetro de escala que será substituído por
δ
√
λ
i
onde este for requerido.
Desta forma, o algoritmo pode ser implementado por meio de iterações
sucessivas a partir de (66), (15), (53), (68), (70), (69) e da distribuição Dirichlet
(60).
4.2 Resultados do ajuste para os três algoritmos
4.2.1 Algoritmo AC (t e normal)
Os resultados obtidos, conforme Figuras 2 e 3, mostram que as estimativas
para os parâmetros de limiar, principalmente, possuem uma autocorrelação muito
forte, dificultando o processo de convergência tanto com função de ligação normal
quanto t de Student acumuladas. Esta análise é coerente com discussões apre-
sentadas por Sörensen et al. (1995), Sörensen & Gianola (2004), Cowles (1996)
e Nandram & Chen (1996), pois, estes autores destacaram o mesmo problema
quando utilizaram o algoritmo AC.
Para os contrastes entre tratamentos com o uso desse algoritmo foi verifi-
cado que estes convergem rapidamente (Figura 11).
A representatividade das categorias quanto à quantidade de elementos tem
relação direta com o processo de convergência. A redução do número de categorias
de nove para cinco, conforme sugestão de Kizilkaya et al. (2003), aumentou a
eficiência do processo de amostragem quanto à convergência. As análises com
nove categorias não são aqui apresentadas.
As estimativas de médias de tratamentos, também, apresentaram forte au-
tocorrelação entre seus valores (Figura 2), indicando a necessidade de uma cadeia
longa para a obtenção de amostras adequadas para fins de inferências.
44
Normal
iterações
ββ1
0 2000 4000 6000 8000 10000
0.5 2.0 3.5
t de Student
iterações
ββ1
0 2000 4000 6000 8000 10000
0.5 2.0 3.5
Normal
iterações
ββ2
0 2000 4000 6000 8000 10000
0.5 2.0 3.5
t de Student
iterações
ββ2
0 2000 4000 6000 8000 10000
0.5 2.0 3.5
Normal
iterações
ββ3
0 2000 4000 6000 8000 10000
0.5 2.0 3.5
t de Student
iterações
ββ3
0 2000 4000 6000 8000 10000
0.5 2.0 3.5
FIGURA 2: Cadeias para médias de tratamentos no algoritmo AC e nas duas fun-
ções de ligação normal de t de Student, acumuladas.
O número de iterações descartadas para eliminar o efeito do valor inicial
arbitrário, isto é, o Burn-in (B) e o salto entre duas iterações consecutivas para que
a amostra final obtida fosse considerada não correlacionada, ou seja, o Thinning
(T) foi obtido por meio do critério de Raftery & Lewis (Tabela 6). Os maiores va-
lores de B (descarte) e T (salto) foram, respectivamente, 248 e 31 para o parâmetro
γ
2
, no modelo com função de ligação normal acumulada e 510 e 24 para γ
3
e σ
2
u
45
Normal
iterações
γγ2
0 2000 4000 6000 8000 10000
0.2 0.8 1.4
t de Student
iterações
γγ2
0 2000 4000 6000 8000 10000
0.5 1.0 1.5
Normal
iterações
γγ3
0 2000 4000 6000 8000 10000
0.5 1.5
t de Student
iterações
γγ3
0 2000 4000 6000 8000 10000
0.5 1.5
Normal
iterações
γγ4
0 2000 4000 6000 8000 10000
1.0 2.0 3.0
t de Student
iterações
γγ4
0 2000 4000 6000 8000 10000
1.5 3.0
FIGURA 3: Cadeias para os parâmetros de limiar no algoritmo AC e modelos nor-
mal e t de Student.
no modelo com função de ligação t de Student acumulada.
Foram realizadas análises para os dados distribuídos em nove categorias
em que os valores de B e T considerados foram, respectivamente, de 5000 e 100
em uma cadeia de 505000 e, portanto, obtendo uma amostra final de 5000 valores
e, ainda assim, não obteve-se convergência. Nesta análise o algoritmo utilizado foi
o AC com distribuição normal para a variável latente.
46
TABELA 6: Valores de Burn-in (B) e Thinning (T) obtidos por meio do diagnós-
tico de Raftery & Lewis, para as médias de tratamentos (β), variância
de blocos (σ
2
u
), variância residuals (δ
2
) e graus de liberdade (ν), nos
algoritmos AC, MC e NC segundo a distribuição da variável latente.
AC MC NC
Normal t Normal t Normal t
Parâmetro B T B T B T B T B T B T
β
1
54 9 48 12 15 5 15 5 3 1 6 2
β
2
45 9 60 10 15 5 6 2 4 2 4 2
β
3
40 8 66 11 9 3 20 5 4 2 3 1
γ
2
45 9 50 10 30 10 128 16 4 2 4 2
γ
3
248 31 117 13 100 20 225 25 2 2 4 2
γ
4
209 19 510 10 85 17 252 28 - - - -
σ
2
u
9 3 96 24 9 3 12 4 6 2 6 2
ν - - 20 5 - - 72 18 - - 78 3
δ
2
- - - - - - - - 6 2 9 3
O maior fator de Brooks, Gelman & Rubin (
R) foi de 1,2584054 também
para γ
2
no modelo normal e de 1,2233068 para γ
3
no modelo t de Student. Se-
gundo Gelman & Rubin (1992), para que haja convergência o valor de
R deve
estar próximo de 1. No entanto, segundo os mesmos autores, um valor próximo de
1,2 pode ser satisfatório para muitos casos. Se
R ≈ 1, então o valor inicial arbi-
trário da cadeia não exerce mais efeito sobre as iterações (Paulino et al., 2003). É
conveniente ressaltar que o valor
R é calculado usando a segunda metade de cada
cadeia.
Em várias outras amostras obtidas esses valores apresentaram outra con-
figuração, isto é, eles possuem oscilação. Porém com um valor de B de alguns
poucos milhares de iterações e um valor de T de algumas poucas dezenas seria o
suficiente para atender às exigências indicadas pelo critério de Raftery & Lewis
em todas as situações.
O fato de considerar para a função G em (2) a distribuição t de Student
47
acumulada, não melhora a convergência do algoritmo em relação ao modelo nor-
mal (Figuras 2 e 3) e os valores de B, T e
R também, não diferem expressivamente
do caso com função de ligação normal acumulada. Neste caso, o modelo é mais
complexo, pois há acréscimo de parâmetros visto que a distribuição t de Student
foi expressa por meio da mistura de uma normal e uma Gama para cada um dos
parâmetros λ
i
, associado às observações. A distribuição Gama para o parâmetro
λ
i
é dada em função do grau de liberdade ν, o qual não tem uma distribuição
“fechada” para a implementação computacional, sendo necessário o uso do algo-
ritmo Metropolis-Hastings para a amostragem de seus valores e isto contribui para
a complexidade do algoritmo.
Os critérios de convergência aplicados e comentados anteriormente e prin-
cipalmente os valores de B e T apontados pelo critério de Raftery & Lewis indicam
a necessidade de cadeias pequenas em contraposição ao que se pode observar por
meio das Figuras 2 e 3, pois verifica-se forte dependência entre as estimativas ob-
tidas em cada iteração. Na Tabela 7 está ilustrada a dependência da média do
tratamento e do parâmetro de limiar que apresentaram maior autocorrelação, tanto
no caso em que a função de ligação foi a distribuição normal acumulada quanto no
caso em que essa função foi a t de Student acumula, conforme Tabela 7. A justifi-
cativa para a apresentação da autocorrelação de apenas dois parâmetros se deve ao
fato de que se a convergência não acontecer para apenas algumas das quantidades
estimadas é suficiente para ilustrar a ineficiência do algoritmo.
Observa-se, ainda, que existe forte autocorrelação para esses parâmetros
principalmente para lag 50 e lag 100 e no modelo t de Student verifica-se forte
autocorrelação também no lag 200. Estes resultados são coerentes com as cadeias
ilustradas nas Figuras 2 e 3.
Conforme valores de ν justifica-se o uso da distribuição t de Student neste
48
TABELA 7: Autocorrelação entre as estimativas para β
1
e para γ
3
no caso em que
a função de ligação foi a distribuição normal acumulada e para β
2
e γ
2
quando a função de ligação foi a distribuição t de Student acumulada,
no algoritmo AC.
Distribuição Parâmetro lag
50 100 200
Normal β
1
0,1761 0,1039 0,0474
γ
3
0,2264 0,1405 0,0637
t β
2
0,3273 0,2451 0,1849
γ
2
0,3420 0,2980 0,2330
iterações
νν
0 2000 4000 6000 8000 10000
0 40 100
νν
Frequência
0 10 20 30 40 50 60
0 4000
FIGURA 4: Cadeia e distribuição das estimativas de ν no algoritmo AC.
algoritmo (Figura 4), pois a maioria das estimativas foram abaixo de 10. No en-
tanto, a dependência entre as estimativas para os parâmetros é maior neste caso.
Este problema é solucionado adotando um grande valor de T.
Estudo de convergência na amostra válida para inferências, isto é, nas ca-
deias sem o efeito do valor inicial arbitrário e estimativas não autocorrelacionadas,
49
foi realizado por meio do critério de Gelman & Rubin. Este estudo, apresen-
tou pequenas variações em torno do valor 1 e, portanto, considera-se que houve
convergência, mesmo neste algoritmo, no qual foram verificadas dificuldades de
convergência a partir do estudo inicial com 10000 iterações. Maiores valores de B
e T dentre todas as cadeias foram, respectivamente, de 6 e 2 quando foi utilizada a
distribuição normal padrão como função de ligação e de 16 e 5 quando se usou a
distribuição t de Student na cadeia referente aos graus de liberdade (ν).
Para a escolha da melhor modelagem foi utilizado o quociente entre os
valores obtidos a partir da expressão (40) com o uso da função de ligação, a dis-
tribuição normal padrão acumulada e da distribuição t de Student acumulada. A
aproximação ao fator de Bayes, dada pela razão mencionada foi de 0,005862994
e tomando o inverso deste valor verifica-se, a partir da Tabela 3, que a evidência é
decididamente favorável ao ajuste por meio da distribuição t de Student.
4.2.2 Algoritmo MC (t e normal)
Observando as Figuras, 2 e 5 e 3 e 6 verifica-se que não houve aumento da
eficiência do algoritmo MC em relação ao AC. Por meio dos valores de B e T indi-
cados para os parâmetros de limiar e apresentados a seguir, verifica-se que houve
um pequeno aumento na eficiênica do algoritmo MC. Os valores de autocorrelação
para estes parâmetros indicam conclusões contrárias pois, em relação ao algoritmo
AC possuem maior dependência entre as estimativas, como é ilustrado para γ
3
nas
Tabelas 7 e 8.
Os maiores valores de B e T, no caso do uso da distribuição normal acumu-
lada em (2) foram, respectivamente, 100 e 20 ambos, para o parâmetro γ
3
. Quando
se usou a distribuição t de Student acumulada em (2), esses valores foram 252 e
28 para o parâmetro γ
4
(Tabela 6).
50
Normal
iterações
ββ1
0 2000 4000 6000 8000 10000
0.5 2.0 3.5
t de Student
iterações
ββ1
0 2000 4000 6000 8000 10000
0.5 2.0 3.5
Normal
iterações
ββ2
0 2000 4000 6000 8000 10000
0.5 2.0 3.5
t de Student
iterações
ββ2
0 2000 4000 6000 8000 10000
0.5 2.0 3.5
Normal
iterações
ββ3
0 2000 4000 6000 8000 10000
0.5 2.0 3.5
t de Student
iterações
ββ3
0 2000 4000 6000 8000 10000
0.5 2.0 3.5
FIGURA 5: Cadeias para médias de tratamentos no algoritmo MC e nas duas fun-
ções de ligação normal de t acumuladas.
Por meio do valor de
R, o algoritmo MC apresenta melhor convergência
em relação ao algoritmo AC com maior valor sendo 1,0685525 para parâmetro γ
2
no modelo normal e de 1,2272103 para os graus de liberdade (ν) no modelo t.
Assim como no algoritmo AC, a contradição entre os valores de B e T
em relação à dependência pode ser verificada por meio dos valores destes e da
autocorrelação indicada na Tabela 8.
51
Normal
iterações
γγ2
0 2000 4000 6000 8000 10000
0.0 0.6 1.2
t de Student
iterações
γγ2
0 2000 4000 6000 8000 10000
0.0 0.6 1.2
Normal
iterações
γγ3
0 2000 4000 6000 8000 10000
0.5 1.5
t de Student
iterações
γγ3
0 2000 4000 6000 8000 10000
0.5 1.5
Normal
iterações
γγ4
0 2000 4000 6000 8000 10000
2.0 3.0 4.0
t de Student
iterações
γγ4
0 2000 4000 6000 8000 10000
2.0 3.0 4.0
FIGURA 6: Cadeias para os parâmetros de limiar no algoritmo MC e modelos
normal e t de Student.
Neste algoritmo também a redução de nove para cinco categorias melho-
rou a convergência nas situações em que foram usadas as distribuições normal e t
acumuladas como funções de ligação.
Em relação ao melhor ajuste de modelos, com o uso das distribuições acu-
muladas normal e t de Student como funções de ligação, em (2), Kizilkaya et al.
(2003) decidiram favoravelmente ao uso da segunda.
52
TABELA 8: Autocorrelação entre as estimativas de β
1
e de γ
4
no caso em que a
função de ligação foi a distribuição normal acumulada e de β
1
e γ
3
quando foi usada a distribuição t de Student acumulada como função
de ligação, no algoritmo MC.
Distribuição Parâmetro lag
50 100 200
Normal β
1
0,3643 0,3276 0,0485
γ
4
0,5453 0,4410 0,0596
t β
1
0,5344 0,5009 0,4279
γ
3
0,8777 0,8042 0,6790
Sabe-se que quando ν → ∞ a distribuição t de Student aproxima-se da
distribuição normal padrão e, assim, se uma estimativa para ν é grande, não se
justifica o uso da distribuição t acumulada em (2). Porém, neste algoritmo, as esti-
mativas para ν, em sua grande maioria, foram estimadas abaixo de 20 o que pode
ser verificado na Figura 7, mas para inferir sobre este parâmetro é necessário uma
cadeia longa, haja vista a necessidade de considerar um intervalo grande entre duas
iterações sucessivas para a obtenção de amostra final não correlacionada (Figura
7).
Neste algoritmo, os valores de
R calculados a partir da amostra válida
para inferências, em todas as cadeias, ficaram dentro do intervalo aceitável para
considerar convergência. Os maiores valores de B e T, também calculdados na
amostra final, foram, respectivamente, de 6 e 2 quando se usou a função de ligação
normal padrão acumulada e de 78 e 3 para a cadeia relativa ao parâmetro ν quando
se utilizou a distribuição t de Student (Tabela 6). Para os demais parâmetros, os
valores de B e T foram similares nas duas modelagens. Portanto, por meio dos
dois critérios verifica-se que houve convergência do algoritmo.
A contrastação dos dois modelos concorrentes, por meio da aproximação
do fator de Bayes, indica que a decisão deve ser tomada com evidência substancial
53
iterações
νν
0 2000 4000 6000 8000 10000
0 50 150
νν
Frequência
0 20 40 60 80 100
0 3000 7000
FIGURA 7: Distribuição das estimativas de ν no algoritmo MC.
(Tabela 3) a favor do uso da distribuição t de Student como função de ligação, já
que a razão entre as médias harmônicas calculadas utilizando (40), no caso em
que a função de ligação foi a normal padrão acumulada e no caso em que foi a
distribuição t de Student, apresentou valor de 0,1862367.
4.2.3 Algoritmo NC (t e normal)
Em contraposição aos algoritmos propostos por Albert & Chib (1993) e por
Cowles (1996) a reparametrização proposta por Nandram & Chen (1996) acelera
expressivamente a convergência, diminuindo a autocorrelação e a necessidade do
descarte de um grande número de valores iniciais para eliminar o efeito do valor
inicial arbitrário.
Sob o aspecto da aceleração da convergência é indiferente o uso da dis-
54
tribuições acumuladas normal ou t de Student em (2), no algoritmo NC também.
A principal diferença entre as duas abordagens é o acréscimo da quantidade de
parâmetros a serem estimados quando se usa a distribuição t e como conseqüência
tem-se maior complexidade na implementação do processo de amostragem.
Com este algoritmo é indiferente, também, a redução do número de cate-
gorias de nove para cinco. Exceto para o parâmetro ν, a convergência do algoritmo
é rápida para todos os parâmetros.
A reparametrização NC consiste na multiplicação dos parâmetros e da va-
riável latente por δ e, desta maneira, para retornar à escala original, basta dividi-los
por este valor.
As Figuras 8 e 9 e os critérios adotados para o estudo de convergência, e
comentados nos dois parágrafos seguintes, ilustram a vantagem do uso do algo-
ritmo NC no modelo normal.
No caso do modelo normal, a indicação das maiores quantidades de B e
T foram, respectivamente, 18 e 6 para o efeito do bloco 35 e de 150 e 25 para ν
quando se utilizou a distribuição t de Student como função de ligação.
Segundo o critério de Gelman & Rubin, as cadeias para todos os parâme-
tros atingiram convergência, dado que o fator determinado por estes autores deve
estar suficientemente próximo de 1 e isto ocorreu nesta análise exceto para o parâ-
metro ν na modelagem com função de ligação t de Student para G em (2).
Esses valores de descarte inicial (B) e pulo (T) são insignificantes quando
comparados à valores que são geralmente adotados. Neste caso, mesmo adotando
um comportamento conservador, e considerando quantidades maiores de valores
para o descarte inicial e para o pulo, não teria uma cadeia tão grande como a
adotada por Sörensen et al. (1995) com 620000 iterações, descarte inicial de 20000
e um valor tomado para a amostra válida a cada 20.
55
Normal
iterações
ββ1
0 2000 4000 6000 8000 10000
0.5 1.5 2.5
t de Student
iterações
ββ1
0 2000 4000 6000 8000 10000
0.5 1.5 2.5
Normal
iterações
ββ2
0 2000 4000 6000 8000 10000
0.0 1.0 2.0
t de Student
iterações
ββ2
0 2000 4000 6000 8000 10000
0.0 1.0 2.0
Normal
iterações
ββ3
0 2000 4000 6000 8000 10000
0.0 1.0 2.0
t de Student
iterações
ββ3
0 2000 4000 6000 8000 10000
0.0 1.0 2.0
FIGURA 8: Cadeias para médias de tratamentos no algoritmo NC e nas duas fun-
ções de ligação normal de t de Student, acumuladas.
Ressalta-se que apesar dos números adotados por Sörensen et al. (1995)
serem uma referência para comparação, a capacidade computacional nos dias atu-
ais os tornam pequenos e, sendo assim, é útil a observação das cadeias e a cons-
tatação de que, a partir de poucas iterações, as amostras para um determinado
parâmetro já são provindas de uma mesma população com o algoritmo NC, e tanto
56
Normal
iterações
γγ2
0 2000 4000 6000 8000 10000
0.0 0.3 0.6
t de Student
iterações
γγ2
0 2000 4000 6000 8000 10000
0.0 0.3 0.6
Normal
iterações
γγ3
0 2000 4000 6000 8000 10000
0.2 0.6 1.0
t de Student
iterações
γγ3
0 2000 4000 6000 8000 10000
0.2 0.6 1.0
Normal
iterações
γγ4
0 2000 4000 6000 8000 10000
1.0 1.6 2.2
t de Student
iterações
γγ4
0 2000 4000 6000 8000 10000
1.0 1.6 2.2
FIGURA 9: Cadeias para os parâmetros de limiar no algoritmo NC e modelos
normal e t.
no modelo normal quanto no modelo t de Student. Além disto, as estimativas pos-
suem pequena autocorrelação (Figuras 8 e 9). Em outras palavras, as cadeias se
estabilizam rapidamente.
É notório a vantagem do algoritmo NC no que se refere à convergência.
Esta afirmativa é confirmada por meio dos critérios utilizados e apontados aqui,
além do comportamento de tendência e dependência das cadeias que podem ser
57
TABELA 9: Autocorrelação entre as estimativas de β
1
e de γ
2
no caso em que a
função de ligação foi a distribuição normal acumulada e de β
1
e de γ
3
quando foi usada a distribuição t de Student acumulada como função
de ligação, no algoritmo NC.
Distribuição Parâmetro lag
5 10 50
Normal β
1
0,1198 0,0436 0,0003
γ
2
0,0750 0,0197 0,0237
t β
1
0,0962 0,0455 -0,0118
γ
3
0,1274 0,0092 0,0068
observadas nas Figuras mencionadas.
Os valores de autocorrelação para os casos mais críticos de convergência,
dentre as médias de tratamentos e dentre os parâmetros de limiar, são indicados na
Tabela 9, assim como nos outros algoritmos. Porém, neste caso, foram considera-
dos intervalos (lags) menores entre duas iterações para o cálculo da autocorrelação
por considerar dispensáveis valores altos.
Os altos valores de ν observados por meio da Figura 10 implicam que a
distribuição t de Student aproxima satisfatoriamente de uma distribuição normal
padrão, o que permite escolher como função de ligação, neste algoritmo, a distri-
buição normal acumulada.
Porém, o parâmetro referente aos graus de liberdade (ν) apresenta proble-
mas de convergências em todos os algoritmos e como em Kizilkaya et al. (2003),
nesta análise também, houve problemas com relação à “mistura”, ou seja, não
houve uma boa aleatorização dos valores amostrados (Figura 10), sendo necessário
um intervalo grande entre duas iterações consecutivas amostradas para a obtenção
de uma amostra final válida para a realização de inferências.
Outros parâmetros, para o algoritmo NC apresentaram rapidez de conver-
gência similar aos citados na Tabela 9, indicando a necessidade de pequenos valo-
58
iterações
νν
0 2000 4000 6000 8000 10000
0 150 300
νν
Frequência
0 50 100 150 200 250 300 350
0 1500 3500
FIGURA 10: Cadeia e histograma para ν no algoritmo NC.
res de T e B.
Na amostra livre do efeito do valor inicial arbitrário e da autocorrelação
encontrou-se valores de
R dentro do aceitável para considerar a convergência do
algoritmo, exceto para a cadeia relativa à ν que foi de 1,3888694. No caso em
que foi utilizada a distribuição normal padrão como função de ligação, o maior e
menor valor de
R foram de 1,00167 e de 0,99980, ambos para efeito de blocos.
Quando se usou a distribuição t de Student, o maior e menor valor de
R foram, res-
59
pectivamente, de 1,00238 e de 0,99983 para média do tratamento 3 (β
3
) e efeito de
bloco. Os maiores valores de B e T foram, respectivamente, 3 e 1 com a primeira
função de ligação e 4 e 2 com a segunda.
Os valores estimados para os graus de liberdade são altos e indicam que a
função de ligação pode ser a distribuição normal padrão acumulada e isto é con-
firmado por meio da aproximação do Fator de Bayes, calculada por meio de da
expressão (40). Este critério permite concluir que a modelagem com o uso da
distribuição normal padrão como função de ligação é ligeiramente superior à t
de Student com o valor da aproximação do fator de Bayes de 1,644553. Porém,
segundo os critérios apresentados na Tabela 3, a evidência de que um modelo é
melhor do que o outro, neste caso, é duvidosa.
Como destacado por Kizilkaya et al. (2003), inferências sobre ν não são
confiáveis dada a dificuldade de convergência do algoritmo para este parâmetro.
Neste trabalho, mesmo após a consideração do Burn-in e do Thinning, verificou-
se que a cadeia para este parâmetro não é estável, com ocorrência de picos. Uma
alternativa para amenizar este problema é o truncamento da distribuição deste pa-
râmetro em 30, pois a distribuição t de Student com graus de liberdade acima deste
valor aproxima-se da distribuição normal padrão.
No algoritmo NC é indiferente o uso da distribuição normal ou t de Student
para a variável latente, o que foi observado a partir do fator de Bayes para os
modelos com essas duas distribuições. Neste algoritmo, a distribuição geradora de
candidatos é a distribuição Dirichlet, a qual, por sua característica, faz com que os
parâmetros threshold fiquem limitados no intervalo [0, 1] e isto modifica a escala
da variável latente, fazendo com que os dois modelos tenham caudas semelhantes.
60
4.2.4 Análise dos contrastes de interesse
Os resultados de contrastes, bem como da correlação intraclasse (ρ) são
apresentados a seguir. Nos algoritmos AC e MC a variância residual foi conside-
rada fixa e igual a 1 (δ
2
= 1).
Em caso de similaridade dessas estatísticas pode se decidir pela equivalên-
cia dos métodos quanto à inferências. Porém, ressalta-se que inferências confiáveis
devem ser realizadas a partir de cadeias estáveis, ou seja, daquelas em que houve
convergência e, por conseqüência, amostras da distribuição condicional conjunta
a posteriori representam uma boa aproximação das distribuições marginais dos
parâmetros.
Na Tabela 10 estão listados os HPDs em todas as situações estudadas para
o contraste associado à regressão linear (C
1
θ), associado ao desvio de regressão
(C
2
θ) e à correlação intraclasse (ρ).
É possível verificar por meio da Tabela 10 que os resultados são similares
em todos os algoritmos, independente da função de ligação utilizada. No estudo
inicial, com 10000 iterações, apenas o algoritmo NC apresentou estabilidade em
suas cadeias, porém com os valores de B e T adotados para a obtenção da amostra
válida para inferências obteve-se convergência em todos os casos. Os HPDs e as
médias apresentadas na Tabela 10 foram calculados a partir da amostra válida.
Apesar da dificuldade de convergência dos algoritmos AC e MC e, prin-
cipalmente para o primeiro, a convergência dos contrastes é evidente. Ilustração
desta situação é apresentada na Figura 11. Os valores de B e T, iguas a 6 e 2,
respectivamente, foram suficientes para eliminar o efeito do valor inicial arbitrário
e para a obteção de uma amostra final não autocorrelacionada. A maior autocor-
relação foi detectada para o efeito da regressão linear com função de ligação t de
Student acumulada com 0,4411 no lag 1 e 0,0894 no lag 5, indicando que em todos
61
TABELA 10: HPDs para funções dos parâmetros associadas à regressão linear
(C
1
β), ao desvio de regressão (C
2
β) e à correlação intraclasse (ρ).
Algoritmo Distribuição Contraste HPD Média
AC Normal C
1
β [-0,5208; 0,0212] -0,2581
t C
1
β [-0,6716; 0,0534] -0,2981
Normal C
2
β [-0,8134; 0,1346] -0,3355
t C
2
β [-1,0329; 0,2129] -0,4190
Normal ρ [ 0,4515; 0,7632] 0,6029
t ρ [ 0,5195; 0,8522] 0,6921
MC Normal C
1
β [-0,5473; 0,0156] -0,2562
t C
1
β [-0,6470; 0,0454] -0,2887
Normal C
2
β [-0,8214; 0,1389] -0,3345
t C
2
β [-1,0243; 0,2252] -0,4035
Normal ρ [ 0,4447; 0,7458] 0,5997
t ρ [ 0,5045; 0,8416] 0,6788
NC Normal C
1
β [-0,6470; 0,0454] -0,2887
t C
1
β [-0,2615; 0,0851] -0,0866
Normal C
2
β [-1,0243; 0,2252] -0,4034
t C
2
β [-0,4212; 0,1829] -0,1208
Normal ρ [ 0,5168; 0,7755] 0,6439
t ρ [ 0,5100; 0,7698] 0,6425
os contrastes a convergência é obtida com um número pequeno de iterações.
4.3 Comparações com a ANAVA e considerações gerais
A análise de variância usual, para os dados sem transformação e com trans-
formação logarítmica, indica que a porcentagem de sacarose influencia significa-
tivamente a nota atribuida à cor da banana da terra desidratada (Tabela 11) o que
não ocorreu quando foi aplicada a transformação de Box-Cox. Foi detectado efeito
linear significativo por meio do teste F, a 5% de probabilidade, apenas no caso da
transformação logarítmica.
O poder da análise reduziu, no caso em que foi usada a transformação
de Box-Cox, pois, nesta análise não foi detectado influência significativa dos trata-
62
Normal
iterações
C1ββ
0 2000 6000 10000
−0.2 0.2 0.6
t de Student
iterações
C1ββ
0 2000 6000 10000
−0.2 0.2 0.6
Normal
iterações
C2ββ
0 2000 6000 10000
−1.0 −0.5 0.0 0.5
t de Student
iterações
C2ββ
0 2000 6000 10000
−1.5 −0.5 0.5
FIGURA 11: Cadeias para os contrastes associados à regressão linear e ao desvio
de regressão no algoritmo AC e funções de ligação normal e t de
Student, acumuladas.
mentos na variável cor (Tabela 11). O valor de λ para a transformação de Box-Cox
foi de 3,5.
Em todos os três casos, conforme valor p da Tabela 11, a falta de ajuste
foi não significativo. Isto é coerente com o contraste C
2
β estimado a partir das
cadeias para a média de cada tratamento pois, o HPD para este contraste não inclui
o valor zero e pode-se concluir então que o desvio de regressão é não significativo.
63
TABELA 11: Quadro-resumo da análise da variância para a variável cor sem
transformação, com transformação logarítmica e com transforma-
ção Box-Cox.
Tipo de transformação
Sem Logarítmica Box-Cox
Fonte Variação GL F valor p F valor p F valor p
Bloco 35 3,18 0,0000 3,19 0,0000 3,18 0,0000
Tratamento 2 3,19 0,0471 3,80 0,0272 2,38 0,0999
Efeito linear 1 3,58 0,0626 4,06 0,0476 2,86 0,0954
Desvio de regressão 1 2,81 0,0984 3,53 0,0644 1,90 0,1721
A hipótese de normalidade dos resíduos, segundo o teste de Shapiro-Wilk
não foi aceita em nenhum dos casos. As probabilidades de significância foram
respectivamente 0,0004421, 0,000007276 e 0,001082 nas análises com os dados
sem transformação, com transformação logarítmica e com transformação de Box-
Cox. A ausência de normalidade dos resíduos pode ser observada também, na
Figura 12.
Em resumo, nenhuma das análise de variâncias pode ser considerada vá-
lida, mas a mais próxima de ser válida (usando a transformação Box-Cox) apre-
senta conclusões semelhantes às dos algoritmos propostos.
64
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
● ●
●
●
●
●
●
●● ●
●●
●
●●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●●●● ●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●●
●
●
●●
●●
●
●●●
●
●
●●
−2 −1 0 1 2
−2 −1 0 1 2
Normal Q−Q Plot (a)
Quantis Teóricos
Resíduos
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●●
●
●●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●●
●●
●
● ●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●●
●
● ●●
●
●
●
●
−2 −1 0 1 2
−0.5 0.0 0.5
Normal Q−Q Plot (b)
Quantis Teóricos
Resíduos
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●●
●
●●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●●
●
● ●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●●
●
● ●●
●
●
●
●
−2 −1 0 1 2
−6 −4 −2 0 2 4 6
Normal Q−Q Plot (c)
Quantis Teóricos
Resíduos
FIGURA 12: Gráfico de resíduos a partir do modelo ajustado para a variável nota
atribuída à cor da banana da terra desidratada, sem transformação
(a), com transformação logarítmica (b) e com transformação de Box-
Cox (c).
65
5 CONCLUSÕES
• As adaptações para modelos mistos do algoritmo NC com funções de liga-
ção normal padrão e t de Student são as mais eficientes e podem ser reco-
mendadas.
• A distribuição t de Student é em geral uma função de ligação mais flexível,
mesmo que no exemplo para o algoritmo NC seja indiferente o seu uso em
alternativa à distribuição normal.
• A análise bayesiana de parâmetros de limiar é recomendável na modelagem
de dados discretos ordinais provenientes da análise sensorial de alimentos.
66
REFERÊNCIAS BIBLIOGRÁFICAS
ALBERT,J. H.; CHIB, S. Bayesian analysis of binary and polychotomous
response data. Journal of the American Statistical Association, Washington, v.
88, n. 422, p. 669-679, June 1993.
BOX, G.E.P; COX, D.R. An analysis of transformations (with discussion).
Journal of the Royal Statistical Society, Series B, London, v. 26, n. 2, p.
211-252, Apr. 1964.
CHIB, S. Markov chain Monte Carlo methods: Computation and Inference. In:
HECKMAN, J. J. and LEAMER E. (Eds.). Handbook of Econometrics:
Volume 5. North-Holland: Elsevier Science, 2001. p. 3569-3649.
COWLES, M. K. Accelerating Monte Carlo Markov chain convergence for
cumulative link generalized linear models, Statistics and Computing,
Netherlands: Springer, v.6, n. 2 p. 101-111, June 1996.
DUTCOSKY, S. D. Análise sensorial de alimentos, Curitiba:Champagnat, 1996.
123p.
GELFAND, A.E.; SMITH A.F.M. Sampling-based approaches to calculating
marginal densities, Journal of the American Statistical Association,
Washington, 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. and RUBIN, D. B. Inference from iterative simulation using
multiple sequences. Statistical Science, Hayward, v. 7, n. 4, p. 457-72, Nov.
1992.
HASTINGS, W. K. Monte Carlo sampling methods using Markov chains and
their applications. Biometrika, London, v. 57, n. 1, p. 97-109, Apr. 1970.
JEFFREYS, H. Theory of probability. Oxford: Claredon Press, 1961. 470 p.
KEMPTHORNE, O.; An introduction to genetic statistic. New York: John
Wiley and Sons Inc., 1966. 545 p.
67
KIZILKAYA, K.; BANKS, B.D.; CARNIER, P.; ALBERA, A.; BITTANTE, G.;
TEMPELMAN, R.J. Bayesian inference strategies for the prediction of genetic
merit using threshold models with an application to calving ease scores in Italian
Piemontese cattle. Journal of Animal Breeding and Genetics, Jokioinen, v.
119, n. 4, p. 209-220, Aug. 2002.
KIZILKAYA, K.; CARNIER, P.; ALBERA, A.; BITTANTE, G.; TEMPELMAN,
R. J. Cumulative t-link threshold models for genetic analysis of calving ease
scores, Genetics Selection Evolution , Les Ulis, v. 35, p. 489-512, Mar. 2003.
LIU, J. S.; WONG, W. H.; KONG, A. Covariance structure of the Gibbs sampler
with applications to the comparison of estimators and augmentation schemes,
Biometrika, London, v. 81, n. 1, p. 27-40, Mar. 1994.
McCULLAGH, P. Regression models for ordinal data, Journal of the Royal
Statistical Society, Oxford, v. 42, n. 2, p. 109-142, Feb. 1980.
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, p. 1087-1092, June 1953.
NANDRAM, B.; CHEN, M. H. Reparameterizing the generalized linear model to
accelerate gibbs sampler convergence, Journal of Statistical Computation and
Simulation, Abingdon, v. 54, p.129-144, Apr. 1996.
PAULINO, C. D.; TURKMAN, M. A. A.; MURTEIRA, B. Estatística
Bayesiana, Lisboa: Fundação Calouste Gulbenkian, 2003, 446 p.
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, p. 493-497. 1992.
R DEVELOPMENT CORE TEM. R: A language and environment for statistical
computing. R Foundation for Statistical Computing, Vienna, Austria. Disponível
em: <http://www.R-project.org.>. Acesso em: 02 de Jul. 2007.
SHAPIRO, S. S.; WILK, M. B. An analysis of variance test for normality
(complete samples), Biometrika, London, v. 52, n. 3 and 4, p.591-611, 1965.
68
SÖRENSEN, D.A.; ANDERSEN S.; GIANOLA D.; KORSGAARD I. Bayesian
inference in threshold models using Gibbs sampling, Genetics Selection
Evolution, Les Ulis, v. 27, n. 3, p. 229-249, 1995.
SÖRENSEN, D.; GIANOLA, D. Likelihood, bayesian and MCMC methods in
quantitative genetics. United States of America: Springer, 2004. 740 p.
WRIGHT, S.; An analysis of variability in number of digits in an inbred strain of
guinea pigs, Genetics, Pittsburgh, v. 19, p. 506-536, Nov. 1934.
69
ANEXO
PROGRAMAS Páginas
PROGRAMA 1 Rotina R para implementação do algoritmo AC no caso em
que a variavel latente tem distribuição normal. . . . . . . . . . . . 71
PROGRAMA 2 Rotina R para implementação do algoritmo AC no caso em
que a variavel latente tem distribuição t de Student. . . . . . . . 75
PROGRAMA 3 Rotina R para implementação do algoritmo MC no caso em
que a variavel latente tem distribuição normal. . . . . . . . . . . . 79
PROGRAMA 4 Rotina R para implementação do algoritmo MC no caso em
que a variavel latente tem distribuição t de Student. . . . . . . . 84
PROGRAMA 5 Rotina R para implementação do algoritmo NC no caso em
que a variavel latente tem distribuição normal. . . . . . . . . . . . 90
PROGRAMA 6 Rotina R para implementação do algoritmo NC no caso em
que a variavel latente tem distribuição t de Student. . . . . . . . 95
70
Programa 1. Rotina R para implementação do algoritmo AC no caso em que a
variavel latente tem distribuição normal.
################################################
## Algoritmo AC #
## Variável latente com distribuição normal #
################################################
dados <- read.table("dados.txt",h=T)
attach(dados)
###### Definição do tamanho da cadeia de saida #
amostra <- 5000
burn <- 5000
jump <- 50
iter <- jump
*
amostra+burn
################################################
library(VGAM)
library(MASS)
library(mvtnorm)
y <- cor
n <- length(y)
#### reorganização dos dados em 5 categorias ##
for(i in 1:n)
{
if (y[i]<=5) y[i]<-1
}
for(i in 1:n)
{
if (y[i] == 6) y[i]<- 2
}
for(i in 1:n)
{
if (y[i] == 7) y[i]<- 3
}
for(i in 1:n)
{
if (y[i] == 8) y[i]<- 4
}
for(i in 1:n)
{
if (y[i] == 9) y[i]<- 5
}
## Matrizes de incidência de tratamentos e blocos #
trat <- c(rep(1,36),rep(2,36),rep(3,36))
bloco <- dados$bloco
XT <-matrix(0,n,3)
71
for(i in 1:n)
{
XT[i,trat[i]]<-1
}
XB <-matrix(0,n,35)
for(i in 1:n)
{
ifelse(bloco[i]==36, XB[i,] <- -1, XB[i,bloco[i]] <- 1)
}
W <- cbind(XT,XB)
######## Valores iniciais arbitrários ##############
L <- c(1:n)
for(i in 1:n)
{
L[i]<- qnorm(pnorm(0)+pnorm(0)
*
runif(1))
}
theta <- ginv(t(W)%
*
%W)%
*
%t(W)%
*
%L
e <- L - W%
*
%theta
A <- diag(35) ## constante
ve <- 1 ## constante
vu <- 3
taunovo <- c(0, 0.2, 0.7, 1.1)
tau <- taunovo
cont <- 0
delta2 <- 1
#### objetos necessários para o cálculo da verossimilhança
vero <- rep(1,n)
verossim <- NULL
########## prioris ###############################
se <- 1000
ru <- 3
su <- 5
################## Amostrador de Gibbs ###########
for (i in 1:iter)
{
############## condicional para theta ############
V <- rbind(cbind(se
*
diag(3) , matrix(0,3,35)),
cbind(matrix(0,35,3), vu
*
A))
S <- ginv(V)
theta <- ginv(t(W)%
*
%W + S)%
*
%t(W)%
*
%L
var <- ve
*
ginv(t(W)%
*
%W+ S)
thetaest <- rmvnorm(1,mean = theta,sigma = var, method = "chol")
theta <- t(thetaest)
beta <- theta[1:3]
u <- theta[4:38]
72
##### condicional para variancia de blocos ########
c1 <- (35+2
*
ru)/2
c2 <- (t(u)%
*
%ginv(A)%
*
%u + 2
*
su)/2
vu <- rgamma(1,c1,c2)
vu <- as.real(1/vu)
############ condicional para L ##################
m <- W%
*
%theta
for(k in 1:n)
{
if(y[k]==1)
{
L[k] <- qnorm(pnorm(tau[1],m[k],ve)
*
runif(1),mean=m[k],sd=sqrt(ve))
vero[k] <- pnorm((tau[1] - m[k])/sqrt(ve))
}
if(y[k]==2)
{
L[k] <- qnorm(pnorm(tau[1],m[k],ve)+(pnorm(tau[2],m[k],ve)-
pnorm(tau[1],m[k],ve))
*
runif(1),mean=m[k],sd=sqrt(ve))
vero[k] <- pnorm((tau[2] - m[k])/sqrt(ve)) -
pnorm((tau[1] - m[k])/sqrt(ve))
}
if(y[k]==3)
{
L[k] <- qnorm(pnorm(tau[2],m[k],ve)+(pnorm(tau[3],m[k],ve)-
pnorm(tau[2],m[k],ve))
*
runif(1),mean=m[k],sd=sqrt(ve))
vero[k] <- pnorm((tau[3] - m[k])/sqrt(ve)) -
pnorm((tau[2] - m[k])/sqrt(ve))
}
if(y[k]==4)
{
L[k] <- qnorm(pnorm(tau[3],m[k],ve)+(pnorm(tau[4],m[k],ve)-
pnorm(tau[3],m[k],ve))
*
runif(1),mean=m[k],sd=sqrt(ve))
vero[k] <- pnorm((tau[4] - m[k])/sqrt(ve)) -
pnorm((tau[3] - m[k])/sqrt(ve))
}
if(y[k]==5)
{
L[k] <- qnorm(pnorm(tau[4],m[k],ve) +
(1-pnorm(tau[4],m[k],ve))
*
runif(1),mean=m[k],sd=sqrt(ve))
vero[k] <- (1 - pnorm( - m[k]/sqrt(ve)))
}
}
############# Cálculo da verossimilhança #######
verossim <- prod(vero)
73
####### condicional para os par. de limiar ####
min <- as.vector(tapply(L,y,min))
max <- as.vector(tapply(L,y,max))
for (j in 2:4)
{
tau[j] <-runif(1,(max[j]),(min[j+1]))
}
###### Gera a saída no arquivo cadeia.txt ######
saida <- cbind(t(theta),t(tau),vu,verossim)
colunas <- length(saída)
if(i > burn && (i-burn)%%jump == 0)
{
write(saida, file="cadeia.txt", ncol=colunas, append=TRUE)
}
}
74
Programa 2. Rotina R para implementação do algoritmo AC no caso em que a
variavel latente tem distribuição t de Student.
################################################
## Algoritmo AC #
## Variável latente com distribuição t #
################################################
dados <-read.table("dados.txt",h=T)
attach(dados)
###### Definição do tamanho da cadeia de saída #
amostra <- 5000
burn <- 5000
jump <- 50
iter <- jump
*
amostra+burn
################################################
library(VGAM)
library(MASS)
library(mvtnorm)
y <- dados$cor
n <- length(y)
#### reorganização dos dados em 5 categorias ##
for(i in 1:n)
{
if (y[i]<=5) y[i]<-1
}
for(i in 1:n)
{
if (y[i] == 6) y[i]<- 2
}
for(i in 1:n)
{
if (y[i] == 7) y[i]<- 3
}
for(i in 1:n)
{
if (y[i] == 8) y[i]<- 4
}
for(i in 1:n)
{
if (y[i] == 9) y[i]<- 5
}
## Matrizes de incidência de tratamentos e blocos #
trat <- c(rep(1,36),rep(2,36),rep(3,36))
bloco <- dados$bloco
XT <-matrix(0,108,3)
for(i in 1:108)
75
{
XT[i,trat[i]]<-1
}
XB <-matrix(0,108,35)
for(i in 1:108)
{
ifelse(bloco[i]==36, XB[i,] <- -1, XB[i,bloco[i]] <- 1)
}
W <- cbind(XT,XB)
######## Valores iniciais arbitrários ##############
L <- c(1:n)
for(i in 1:n)
{
L[i]<- qnorm(pnorm(0)+pnorm(0)
*
runif(1))
}
theta <- ginv(t(W)%
*
%W)%
*
%t(W)%
*
%L
e <- L - W%
*
%theta
A <- diag(35) ## constante
nu <- 10
vu <- 0.5
ve <- 1 ## constante
taunovo <- c(0, 0.4, 1.0, 1.4)
tau <- taunovo
pnuc <- 0.5
pnu <- 0.5
lambda <- c(rep(1,n))
R <- diag(n)
#### objetos necessários para o cálculo da verossimilhança
vero <- rep(1,n)
verossim <- NULL
################### prioris ##################################
se <- 1000
ru <- 3
su <- 5
################ Amostrador Gibbs ############################
for (i in 1:iter)
{
################# condicional para theta ###################
V <- rbind(cbind(se
*
diag(3) , matrix(0,3,35)),
cbind(matrix(0,35,3), vu
*
A))
S <- ve
*
ginv(V)
theta <- ginv(t(W)%
*
%R%
*
%W + S)%
*
%t(W)%
*
%R%
*
%L
var <- ve
*
ginv(t(W)%
*
%R%
*
%W+ S)
thetaest <- rmvnorm(1,mean = theta,sigma = var, method = "chol")
theta <- t(thetaest)
beta <- theta[1:3]
76
u <- theta[4:38]
######### condicional para variancia de blocos ################
c1 <- (35+2
*
ru)/2
c2 <- (t(u)%
*
%ginv(A)%
*
%u + 2
*
su)/2
vu <- rgamma(1,c1,c2)
vu <- as.real(1/vu)
################### Amostra lambda (GS elemento a elemento) ###
for (j in 1:n)
{
lambda[j] <- rgamma(1,(nu+1)/2,(nu+(L[j]-W[j,]%
*
%theta)^2)/2)
}
R <- diag(lambda)
################# Amostra nu (MH) ############################
nuc <- 2
while(nuc < 3)
{
nuc <- rpois(1,nu)
}
pnu <- n
*
((nu/2)
*
log(nu/2)-lgamma(nu/2))+sum(((nu/2)-1)
*
log(lambda) + (-(nu/2)
*
lambda)) - 2
*
log(1+nu)
pnuc <- n
*
((nuc/2)
*
log(nuc/2)-lgamma(nuc/2))+sum(((nuc/2)-1)
*
log(lambda) + (-(nuc/2)
*
lambda)) - 2
*
log(1+nuc)
pnu <- exp(pnu)
pnuc <- exp(pnuc)
dnu <- dpois(nu,nu)
dnuc <- dpois(nuc,nu)
num <- pnuc
*
dnu
denom <- pnu
*
dnuc
ifelse((num > denom),(alfa <- 1),(alfa <- as.real(num/denom)))
teste <- runif(1)
if(teste < alfa)
{
nu <- nuc
}
############## condicional para L ########################
m <- W%
*
%theta
dp <- rep(1,n)
dp <- sqrt(ve/lambda)
for(k in 1:n)
{
if(y[k]==1)
{
77
L[k] <- qnorm(pnorm(tau[1],m[k],sd=dp[k])
*
runif(1),mean=m[k],sd=dp[k])
vero[k] <- pnorm((tau[1] - m[k])/dp[k])
}
if(y[k]==2)
{
L[k] <- qnorm(pnorm(tau[1],m[k],sd=dp[k])+(pnorm(tau[2],m[k],sd=dp[k])
-pnorm(tau[1],m[k],sd=dp[k]))
*
runif(1),mean=m[k],sd=dp[k])
vero[k] <- pnorm((tau[2] - m[k])/dp[k]) - pnorm((tau[1] - m[k])/dp[k])
}
if(y[k]==3)
{
L[k] <- qnorm(pnorm(tau[2],m[k],sd=dp[k])+(pnorm(tau[3],m[k],sd=dp[k])
-pnorm(tau[2],m[k],sd=dp[k]))
*
runif(1),mean=m[k],sd=dp[k])
vero[k] <- pnorm((tau[3] - m[k])/dp[k]) - pnorm((tau[2] - m[k])/dp[k])
}
if(y[k]==4)
{
L[k] <- qnorm(pnorm(tau[3],m[k],sd=dp[k])+(pnorm(tau[4],m[k],sd=dp[k])
-pnorm(tau[3],m[k],sd=dp[k]))
*
runif(1),mean=m[k],sd=dp[k])
vero[k] <- pnorm((tau[4] - m[k])/dp[k]) - pnorm((tau[3] - m[k])/dp[k])
}
if(y[k]==5)
{
L[k] <- qnorm(pnorm(tau[4],m[k],sd=dp[k]) +
(1-pnorm(tau[4],m[k],sd=dp[k]))
*
runif(1),mean=m[k],sd=dp[k])
vero[k] <- (1 - pnorm( - m[k]/dp[k]))
}
}
############# Cálculo da verossimilhança #####################
verossim <- prod(vero)
############# Amostragem dos thresholds ######################
min <- as.vector(tapply(L,y,min))
max <- as.vector(tapply(L,y,max))
for(j in 2:4)
{
tau[j] <-runif(1,(max[j]),(min[j+1]))
}
############ Gera a saída no arquivo cadeia.txt #############
saida <- cbind(t(theta),t(tau),nu,vu,verossim)
colunas <- length(saida)
if(i > burn && (i-burn)%%jump == 0)
{
write(saida, file="cadeia.txt", ncol=colunas, append=TRUE)
}
}
78
Programa 3. Rotina R para implementação do algoritmo MC no caso em que a
variavel latente tem distribuição normal.
################################################
## Algoritmo MC #
## Variável latente com distribuição normal #
################################################
dados <- read.table("dados.txt",h=T)
attach(dados)
###### Definição do tamanho da cadeia de saída #
amostra <- 5000
burn <- 5000
jump <- 50
iter <- jump
*
amostra+burn
################################################
library(VGAM)
library(MASS)
library(mvtnorm)
y <- cor
n <- length(y)
#### reorganização dos dados em 5 categorias ##
for(i in 1:n)
{
if (y[i]<=5) y[i]<-1
}
for(i in 1:n)
{
if (y[i] == 6) y[i]<- 2
}
for(i in 1:n)
{
if (y[i] == 7) y[i]<- 3
}
for(i in 1:n)
{
if (y[i] == 8) y[i]<- 4
}
for(i in 1:n)
{
if (y[i] == 9) y[i]<- 5
}
## Matrizes de incidência de tratamentos e blocos #
trat <- c(rep(1,36),rep(2,36),rep(3,36))
bloco <- dados$bloco
XT <-matrix(0,n,3)
for(i in 1:n)
79
{
XT[i,trat[i]]<-1
}
XB <-matrix(0,n,35)
for(i in 1:n)
{
ifelse(bloco[i]==36, XB[i,] <- -1, XB[i,bloco[i]] <- 1)
}
W <- cbind(XT,XB)
## objetos para o cálculo da verossimilhança ##
vero <- rep(1,n)
verossim <- NULL
######## Valores iniciais arbitrários ##############
ve <- 1
vu <- 3
tau <- c( 0, 0.2, 0.7, 1.0, 1000)
taunovo <- c( 0, 0.1, 0.9, 1.2, 1000)
L <- c(1:n)
for(i in 1:n)
{
L[i]<- qnorm(pnorm(0)+pnorm(0)
*
runif(1))
}
############ Número de categorias #############
K <- 5
########### desvio padrão para a candidata ###
dp <- 0.15
d <- rep(1,iter) # vetor para atualização do dp
############## prioris #########################
se <- 1000
ru <- 3
su <- 5
####### Inicia a rotina do amostrador Gibbs ###
for (i in 1:iter)
{
m <- W%
*
%theta
############ condicional para theta ##########
V <- rbind(cbind(se
*
diag(3) , matrix(0,3,35)),
cbind(matrix(0,35,3), vu
*
A))
S <- ve
*
ginv(V)
theta <- ginv(t(W)%
*
%W + S)%
*
%t(W)%
*
%L
var <- ve
*
ginv(t(W)%
*
%W+ S)
80
thetaest <- rmvnorm(1,mean = theta,sigma = var, method = "chol")
theta <- t(thetaest)
beta <- theta[1:3]
u <- theta[4:38]
######## condicional para variancia de blocos###
c1 <- (35+2
*
ru)/2
c2 <- (t(u)%
*
%ginv(A)%
*
%u + 2
*
su)/2
vu <- rgamma(1,c1,c2)
vu <- as.real(1/vu)
############# Amostragem dos thresholds #########
for(j in 2:(K-1))
{
area1 <- pnorm(taunovo[j-1],mean=tau[j],sd=dp)
area2 <- pnorm(tau[j+1],mean=tau[j],sd=dp)-
pnorm(taunovo[j-1],mean=tau[j],sd=dp)
area <- area1+runif(1)
*
area2
taunovo[j] <- qnorm(area,mean=tau[j],sd=dp)
}
p1<-rep(0,K-2)
p2<-rep(0,K-2)
p3<-rep(0,K-2)
p4<-rep(0,K-2)
for(j in 1:(K-2))
{
p1[j] <- pnorm((tau[j+2] -tau[j+1])/dp)
}
for(j in 1:(K-2))
{
p2[j] <- pnorm((taunovo[j] -tau[j+1])/dp)
}
for(j in 1:(K-2))
{
p3[j] <- pnorm((taunovo[j+2] -taunovo[j+1])/dp)
}
for(j in 1:(K-2))
{
p4[j] <- pnorm((tau[j] -taunovo[j+1])/dp)
}
for(k in 1:n)
{
if(y[k]==1)
{
ftau[k] <- pnorm(tau[1] - m[k])
ftaunovo[k] <- pnorm(taunovo[1] - m[k])
}
if(y[k]==2)
{
ftau[k] <- pnorm(tau[2] - m[k]) - pnorm(tau[1] - m[k])
ftaunovo[k] <- pnorm(taunovo[2] - m[k]) - pnorm(taunovo[1] - m[k])
81
}
if(y[k]==3)
{
ftau[k] <- pnorm(tau[3] - m[k]) - pnorm(tau[2] - m[k])
ftaunovo[k] <- pnorm(taunovo[3] - m[k]) - pnorm(taunovo[2] - m[k])
}
if(y[k]==4)
{
ftau[k] <- pnorm(tau[4] - m[k]) - pnorm(tau[3] - m[k])
ftaunovo[k] <- pnorm(taunovo[4] - m[k]) - pnorm(taunovo[3] - m[k])
}
if(y[k]==5)
{
ftau[k] <- 1 - pnorm(tau[4] - m[k])
ftaunovo[k] <- 1 - pnorm(taunovo[4] - m[k])
}
}
alfa <- exp(sum(log((p1-p2))-log((p3-p4)))+sum(log(ftaunovo)-log(ftau)))
teste <- runif(1)
R <- min(1,alfa)
if(teste < R)
{
tau <- taunovo
############### condicional para L #############
for(k in 1:n)
{
if(y[k]==1)
{
L[k] <-qnorm(pnorm(tau[1],m[k])
*
runif(1),mean=m[k])
vero[k] <-pnorm((tau[1] - m[k])/sqrt(ve))
}
if(y[k]==2)
{
L[k] <-qnorm(pnorm(tau[1],mean=m[k])+(pnorm(tau[2],mean=m[k])
-pnorm(tau[1],mean=m[k]))
*
runif(1),mean=m[k])
vero[k] <-pnorm((tau[2] - m[k])/sqrt(ve)) -
pnorm((tau[1] - m[k])/sqrt(ve))
}
if(y[k]==3)
{
L[k] <-qnorm(pnorm(tau[2],mean=m[k])+(pnorm(tau[3],mean=m[k])
-pnorm(tau[2],mean=m[k]))
*
runif(1),mean=m[k])
vero[k] <-pnorm((tau[3] - m[k])/sqrt(ve)) -
pnorm((tau[2] - m[k])/sqrt(ve))
}
if(y[k]==4)
{
L[k] <-qnorm(pnorm(tau[3],mean=m[k])+(pnorm(tau[4],mean=m[k])
-pnorm(tau[3],mean=m[k]))
*
runif(1),mean=m[k])
vero[k] <-pnorm((tau[4] - m[k])/sqrt(ve)) -
pnorm((tau[3] - m[k])/sqrt(ve))
}
82
if(y[k]==5)
{
L[k] <-qnorm(pnorm(tau[4],mean=m[k]) +
(1-pnorm(tau[4],mean=m[k]))
*
runif(1),mean=m[k])
vero[k] <-(1 - pnorm( - m[k]/sqrt(ve)))
}
}
}
############# Cálculo da verossimilhança #####################
verossim <- prod(vero)
######## atualização do desvio padrão da candidata ##########
d[i] <- c(tau[2])
ifelse(i > 20, dp <- sd(d[(i-20):i]), dp <- 0.25)
################# Gera a saída no arquivo cadeia.txt #########
saida <- cbind(t(theta),t(tau),vu,dp,verossim)
colunas <- length(saida)
if(i > burn && (i-burn)%%jump == 0)
{
write(saida, file="cadeia.txt", ncol=colunas, append=TRUE)
}
}
83
Programa 4. Rotina R para implementação do algoritmo MC no caso em que a
variavel latente tem distribuição t de Student.
################################################
## Algoritmo MC #
## Variável latente com distribuição t #
################################################
dados<-read.table("dados.txt",h=T)
###### Definição do tamanho da cadeia de saída #
amostra <- 5000
burn <- 5000
jump <- 50
iter <- jump
*
amostra+burn
################################################
library(VGAM)
library(MASS)
library(mvtnorm)
attach(dados)
y <- cor
n <- length(y)
#### reorganização dos dados em 5 categorias ##
for(i in 1:n)
{
if (y[i]<=5) y[i]<-1
}
for(i in 1:n)
{
if (y[i] == 6) y[i]<- 2
}
for(i in 1:n)
{
if (y[i] == 7) y[i]<- 3
}
for(i in 1:n)
{
if (y[i] == 8) y[i]<- 4
}
for(i in 1:n)
{
if (y[i] == 9) y[i]<- 5
}
## Matrizes de incidência de tratamentos e blocos #
trat <- c(rep(1,36),rep(2,36),rep(3,36))
XT <- matrix(0,108,3)
for(i in 1:108)
{
84
XT[i,trat[i]]<-1
}
XB <- matrix(0,108,35)
for(i in 1:108)
{
ifelse(bloco[i]==36, XB[i,] <- -1, XB[i,bloco[i]] <- 1)
}
W <- cbind(XT,XB)
######## Valores iniciais arbitrários ##############
theta <- ginv(t(W)%
*
%W)%
*
%t(W)%
*
%L
e <- L - W%
*
%theta
A <- diag(35)
ftau <- rep(1,n)
ftaunovo <- ftau
vu <- 3
ve <- 1
tau <- c( 0, 0.3, 0.9, 1.5, 1000)
taunovo <- c( 0, 0.6, 0.8, 2.0, 1000)
nu <- 10
pnu <- 0.5
pnuc <- 0.5
lambda <- c(rep(1,n))
R <- diag(n)
L <- c(1:n)
for(i in 1:n)
{
L[i]<- qnorm(pnorm(0)+pnorm(0)
*
runif(1))
}
############ Número de categorias #############
K <- 5
########### desvio padrão para a candidata ###
dpc <- 0.15
dc <- rep(1,iter) # vetor para atualização do dpc
############## prioris #########################
se <- 1000
ru <- 3
su <- 5
## objetos para o cálculo da verossimilhança ##
vero <- rep(1,n)
verossim <- NULL
####### Inicia a rotina do amostrador Gibbs ###
for (i in 1:iter)
{
m <- W%
*
%theta
85
####### condicional para beta ####
V <- rbind(cbind(se
*
diag(3) , matrix(0,3,35)),
cbind(matrix(0,35,3), vu
*
A))
S <- ve
*
ginv(V)
theta <- ginv(t(W)%
*
%R%
*
%W + S)%
*
%t(W)%
*
%R%
*
%L
var <- ve
*
ginv(t(W)%
*
%R%
*
%W+ S)
thetaest <- rmvnorm(1,mean = theta,sigma = var, method = "chol")
theta <- t(thetaest)
beta <- theta[1:3]
u <- theta[4:38]
### condicional para variancia de blocos ###
c1 <- (35+2
*
ru)/2
c2 <- (t(u)%
*
%ginv(A)%
*
%u + 2
*
su)/2
vu <- rgamma(1,c1,c2)
vu <- as.real(1/vu)
############# Amostragem dos thresholds #############
for(j in 2:(K-1))
{
area1 <- pnorm(taunovo[j-1],mean=tau[j],sd=dpc)
area2 <- pnorm(tau[j+1],mean=tau[j],sd=dpc)-
pnorm(taunovo[j-1],mean=tau[j],sd=dpc)
area <- area1+runif(1)
*
area2
taunovo[j] <- qnorm(area,mean=tau[j],sd=dpc)
}
p1<-rep(0,K-2)
p2<-rep(0,K-2)
p3<-rep(0,K-2)
p4<-rep(0,K-2)
for(j in 1:(K-2))
{
p1[j] <- pnorm((tau[j+2] -tau[j+1])/dpc)
}
for(j in 1:(K-2))
{
p2[j] <- pnorm((taunovo[j] -tau[j+1])/dpc)
}
for(j in 1:(K-2))
{
p3[j] <- pnorm((taunovo[j+2] -taunovo[j+1])/dpc)
}
for(j in 1:(K-2))
{
p4[j] <- pnorm((tau[j] -taunovo[j+1])/dpc)
}
########## Amostra lambda (GS elemento a elemento) ###
for (j in 1:n)
{
86
lambda[j] <- rgamma(1,(nu+1)/2,(nu+(L[j]-W[j,]%
*
%theta)^2)/2)
}
R <- diag(lambda)
for(k in 1:n)
{
m <- W%
*
%theta
dp <- rep(1,n)
dp <- sqrt(ve/lambda)
if(y[k]==1)
{
ftau[k] <- pnorm((tau[1] - m[k])/dp[k])
ftaunovo[k] <- pnorm((taunovo[1] - m[k])/dp[k])
}
if(y[k]==2)
{
ftau[k] <- pnorm((tau[2] - m[k])/dp[k]) -
pnorm((tau[1] - m[k])/dp[k])
ftaunovo[k] <- pnorm((taunovo[2] - m[k])/dp[k]) -
pnorm((taunovo[1] - m[k])/dp[k])
}
if(y[k]==3)
{
ftau[k] <- pnorm((tau[3] - m[k])/dp[k]) -
pnorm((tau[2] - m[k])/dp[k])
ftaunovo[k] <- pnorm((taunovo[3] - m[k])/dp[k]) -
pnorm((taunovo[2] - m[k])/dp[k])
}
if(y[k]==4)
{
ftau[k] <- pnorm((tau[4] - m[k])/dp[k]) -
pnorm((tau[3] - m[k])/dp[k])
ftaunovo[k] <- pnorm((taunovo[4] - m[k])/dp[k]) -
pnorm((taunovo[3] - m[k])/dp[k])
}
if(y[k]==5)
{
ftau[k] <- 1 - pnorm((tau[4] - m[k])/dp[k])
ftaunovo[k] <- 1 - pnorm((taunovo[4] - m[k])/dp[k])
}
}
alfa <- exp(sum(log((p1-p2))-log((p3-p4)))+
sum(log(ftaunovo)-log(ftau)))
teste <- runif(1)
D <- min(1,alfa)
if(teste < D)
{
tau <- taunovo
################### Amostra nu (MH) ##############
nuc <- 2
while(nuc < 3)
87
{
nuc <- rpois(1,nu)
}
pnu <- n
*
((nu/2)
*
log(nu/2)-lgamma(nu/2))+sum(((nu/2)-1)
*
log(lambda)+
(-(nu/2)
*
lambda)) - 2
*
log(1+nu)
pnuc <- n
*
((nuc/2)
*
log(nuc/2)-lgamma(nuc/2))+sum(((nuc/2)-1)
*
log(lambda)+
(-(nuc/2)
*
lambda)) - 2
*
log(1+nuc)
pnu <- exp(pnu)
pnuc <- exp(pnuc)
dnu <- dpois(nu,nu)
dnuc <- dpois(nuc,nu)
num <- pnuc
*
dnu
denom <- pnu
*
dnuc
ifelse((num > denom),(alfa <- 1),(alfa <- as.real(num/denom)))
teste <- runif(1)
if(teste < alfa)
{
nu <- nuc
contnu <- contnu+1
}
########## condicional para L ####
for(k in 1:n)
{
if(y[k]==1)
{
L[k] <-qnorm(pnorm(tau[1],m[k],sd=dp[k])
*
runif(1),mean=m[k],sd=dp[k])
vero[k] <-pnorm((tau[1] - m[k])/dp[k])
}
if(y[k]==2)
{
L[k] <-qnorm(pnorm(tau[1],m[k],sd=dp[k])+(pnorm(tau[2],m[k],sd=dp[k])
-pnorm(tau[1],m[k],sd=dp[k]))
*
runif(1),mean=m[k],sd=dp[k])
vero[k] <-pnorm((tau[2] - m[k])/dp[k]) - pnorm((tau[1] - m[k])/dp[k])
}
if(y[k]==3)
{
L[k] <-qnorm(pnorm(tau[2],m[k],sd=dp[k])+(pnorm(tau[3],m[k],sd=dp[k])
-pnorm(tau[2],m[k],sd=dp[k]))
*
runif(1),mean=m[k],sd=dp[k])
vero[k] <-pnorm((tau[3] - m[k])/dp[k]) - pnorm((tau[2] - m[k])/dp[k])
}
if(y[k]==4)
{
L[k] <- qnorm(pnorm(tau[3],m[k],sd=dp[k])+(pnorm(tau[4],m[k],sd=dp[k])
-pnorm(tau[3],m[k],sd=dp[k]))
*
runif(1),mean=m[k],sd=dp[k])
vero[k] <-pnorm((tau[4] - m[k])/dp[k]) - pnorm((tau[3] - m[k])/dp[k])
}
if(y[k]==5)
{
L[k] <- qnorm(pnorm(tau[4],m[k],sd=dp[k]) +
(1-pnorm(tau[4],m[k],sd=dp[k]))
*
runif(1),mean=m[k],sd=dp[k])
vero[k] <-(1 - pnorm( - m[k]/dp[k]))
}
}
88
}
verossim <- prod(vero)
######## atualização do desvio padrão da candidata ###########
dc[i] <- c(tau[3])
ifelse(i > 20, dpc <- sd(dc[(i-20):i]), dpc <- 0.25)
######### Gera a saída no arquivo cadeia.txt #########
saida <- cbind(t(theta),t(tau),vu,nu,dpc,verossim)
colunas <- length(saida)
if(i > burn && (i-burn)%%jump == 0)
{
write(saida, file="cadeia.txt", ncol=colunas, append=TRUE)
}
}
89
Programa 5. Rotina R para implementação do algoritmo NC no caso em que a
variavel latente tem distribuição normal.
################################################
## Algoritmo NC #
## Variável latente com distribuição normal #
################################################
dados <- read.table("dados.txt",h=T)
attach(dados)
###### Definição do tamanho da cadeia de saída #
amostra <- 5000
burn <- 5000
jump <- 50
iter <- jump
*
amostra+burn
################################################
library(VGAM)
library(MASS)
library(mvtnorm)
y <- cor
n <- length(y)
#### reorganização dos dados em 5 categorias ##
for(i in 1:n)
{
if (y[i]<=5) y[i]<-1
}
for(i in 1:n)
{
if (y[i] == 6) y[i]<- 2
}
for(i in 1:n)
{
if (y[i] == 7) y[i]<- 3
}
for(i in 1:n)
{
if (y[i] == 8) y[i]<- 4
}
for(i in 1:n)
{
if (y[i] == 9) y[i]<- 5
}
## Matrizes de incidência de tratamentos e blocos #
trat <- c(rep(1,36),rep(2,36),rep(3,36))
bloco <- dados$bloco
XT <-matrix(0,n,3)
for(i in 1:n)
90
{
XT[i,trat[i]]<-1
}
XB <-matrix(0,n,35)
for(i in 1:n)
{
ifelse(bloco[i]==36, XB[i,] <- -1, XB[i,bloco[i]] <- 1)
}
W <- cbind(XT,XB)
##### objetos para o cálculo da verossimilhança ####
vero <- rep(1,n)
verossim <- NULL
##### objetos para a atualização dos par. de limiar #
nk <- tapply(L,y,length)
ftau <- rep(1,n-nk[1]-nk[5])
ftaunovo <- ftau
######## Valores iniciais arbitrários ##############
vu <- 3
pnovo <- c(0.2,0.4,0.4)
p <- pnovo
taunovo <- c(0, 0.2, 0.6, 1)
tau <- taunovo
delta2 <- 0.5
delta <- sqrt(delta2)
L <- c(1:n)
for(i in 1:n)
{
L[i]<- qnorm(pnorm(0)+pnorm(0)
*
runif(1))
}
############## prioris #########################
se <- 1000
dre <- 20
dse <- 5
ru <- 3
su <- 5
##################### Amostrador de Gibbs ###
for (i in 1:iter)
{
####### condicional para theta ############
V <- rbind(cbind(se
*
diag(3) , matrix(0,3,35)),
cbind(matrix(0,35,3), vu
*
A))
S <- delta2
*
ginv(V)
theta <- ginv(t(W)%
*
%W + S)%
*
%t(W)%
*
%L
var <- delta2
*
ginv(t(W)%
*
%W+ S)
thetaest <- rmvnorm(1,mean = theta,sigma = var, method = "chol")
91
theta <- t(thetaest)
beta <- theta[1:3]
u <- theta[4:38]
### condicional para variancia de blocos ###
c1 <- (35+2
*
ru)/2
c2 <- (t(u)%
*
%ginv(A)%
*
%u + 2
*
su)/2
vu <- rgamma(1,c1,c2)
vu <- as.real(1/vu)
## par. de aceleramento da convergência (reparametrização) ###
delta2 <- 1/rgamma(1,(n+38+5+2
*
dre)/2,(t(L-W%
*
%theta)%
*
%(L-W%
*
%theta)
+ 2
*
dse +t(theta)%
*
%ginv(V)%
*
%theta)/2)
delta <- sqrt(delta2)
############# Amostragem dos par. de limiar #############
nkp <- 0.8
*
nk[2:(5-1)]
pnovo <- c(rdiric(1, shape = nkp ))
for(j in 2:(5-1))
{
taunovo[j]<-sum(pnovo[1:(j-1)])
}
lfp <- nkp
*
log(p)
lfpnovo <- nkp
*
log(pnovo)
m <- W%
*
%theta
for(k in 1:n)
{
if(y[k]==1)
{
ftau[k] <- 1
ftaunovo[k] <- 1
}
if(y[k]==2)
{
ftau[k] <- pnorm((tau[2] - m[k])/delta) -
pnorm((tau[1] - m[k])/delta)
ftaunovo[k] <- pnorm((taunovo[2] - m[k])/delta) -
pnorm((taunovo[1] - m[k])/delta)
}
if(y[k]==3)
{
ftau[k] <- pnorm((tau[3] - m[k])/delta) -
pnorm((tau[2] - m[k])/delta)
ftaunovo[k] <- pnorm((taunovo[3] - m[k])/delta) -
pnorm((taunovo[2] - m[k])/delta)
}
if(y[k]==4)
{
ftau[k] <- pnorm((tau[4] - m[k])/delta) -
pnorm((tau[3] - m[k])/delta)
ftaunovo[k] <- pnorm((taunovo[4] - m[k])/delta) -
92
pnorm((taunovo[3] - m[k])/delta)
}
if(y[k]==5)
{
ftau[k] <- 1
ftaunovo[k] <- 1
}
}
lftau <- log(ftau)
lftaunovo <- log(ftaunovo)
alfa <- exp(sum(lftaunovo-lftau)-sum(lfpnovo-lfp))
teste <- runif(1)
R <- min(1,alfa)
if(teste < R)
{
cont <- cont + 1
p <- pnovo
tau <- taunovo
########## condicional para L ####
for(k in 1:n)
{
if(y[k]==1)
{
L[k] <-qnorm(pnorm(tau[1],m[k],sd=delta)
*
runif(1),mean=m[k],sd=delta)
vero[k]<-pnorm((tau[1] - m[k])/delta)
}
if(y[k]==2)
{
L[k] <-qnorm(pnorm(tau[1],m[k],sd=delta)+(pnorm(tau[2],m[k],sd=delta)
-pnorm(tau[1],m[k],sd=delta))
*
runif(1),mean=m[k],sd=delta)
vero[k]<-pnorm((tau[2] - m[k])/delta) - pnorm((tau[1] - m[k])/delta)
}
if(y[k]==3)
{
L[k] <-qnorm(pnorm(tau[2],m[k],sd=delta)+(pnorm(tau[3],m[k],sd=delta)
-pnorm(tau[2],m[k],sd=delta))
*
runif(1),mean=m[k],sd=delta)
vero[k]<-pnorm((tau[3] - m[k])/delta) - pnorm((tau[2] - m[k])/delta)
}
if(y[k]==4)
{
L[k] <-qnorm(pnorm(tau[3],m[k],sd=delta)+(pnorm(tau[4],m[k],sd=delta)
-pnorm(tau[3],m[k],sd=delta))
*
runif(1),mean=m[k],sd=delta)
vero[k]<-pnorm((tau[4] - m[k])/delta) - pnorm((tau[3] - m[k])/delta)
}
if(y[k]==5)
{
L[k] <-qnorm(pnorm(tau[4],m[k],sd=delta)+
(1 - pnorm(tau[4],m[k],sd=delta))
*
runif(1),mean=m[k],sd=delta)
vero[k]<-(1 - pnorm( - m[k]/delta))
}
}
}
93
############# Cálculo da verossimilhança #####################
verossim <- prod(vero)
####### Gera a saída no arquivo cadeia.txt ####################
saida <- cbind(t(theta),t(tau),delta2,vu,verossim)
colunas <- length(saida)
if(i > burn && (i-burn)%%jump == 0)
{
write(saida, file="cadeia.txt", ncol=colunas, append=TRUE)
}
}
94
Programa 6. Rotina R para implementação do algoritmo NC no caso em que a
variavel latente tem distribuição t de Student.
################################################
## Algoritmo NC #
## Variável latente com distribuição t #
################################################
dados <- read.table("dados.txt",h=T)
attach(dados)
###### Definição do tamanho da cadeia de saída #
amostra <- 5000
burn <- 5000
jump <- 50
iter <- jump
*
amostra+burn
################################################
library(VGAM)
library(MASS)
library(mvtnorm)
y <- cor
n <- length(y)
#### reorganização dos dados em 5 categorias ##
for(i in 1:n)
{
if (y[i]<=5) y[i]<-1
}
for(i in 1:n)
{
if (y[i] == 6) y[i]<- 2
}
for(i in 1:n)
{
if (y[i] == 7) y[i]<- 3
}
for(i in 1:n)
{
if (y[i] == 8) y[i]<- 4
}
for(i in 1:n)
{
if (y[i] == 9) y[i]<- 5
}
## Matrizes de incidência de tratamentos e blocos #
trat <- c(rep(1,36),rep(2,36),rep(3,36))
bloco <- dados$bloco
XT <-matrix(0,n,3)
for(i in 1:n)
95
{
XT[i,trat[i]]<-1
}
XB <-matrix(0,n,35)
for(i in 1:n)
{
ifelse(bloco[i]==36, XB[i,] <- -1, XB[i,bloco[i]] <- 1)
}
W <- cbind(XT,XB)
##### objetos para o cálculo da verossimilhança ####
vero <- rep(1,n)
verossim <- NULL
##### objetos para a atualização dos par. de limiar #
nk <- tapply(L,y,length)
ftau <- rep(1,n-nk[1]-nk[5])
ftaunovo <- ftau
######## Valores iniciais arbitrários ##############
vu <- 3
pnovo <- c(0.2,0.4,0.4)
p <- pnovo
taunovo <- c(0, 0.2, 0.6, 1)
tau <- taunovo
delta2 <- 0.5
delta <- sqrt(delta2)
nu <- 8
pnu <- 0.5
pnuc <- 0.5
nuc <- rpois(1,5)
lambda <- c(rep(1,n))
R <- diag(n)
L <- c(1:n)
for(i in 1:n)
{
L[i]<- qnorm(pnorm(0)+pnorm(0)
*
runif(1))
}
############## prioris #########################
se <- 1000
dre <- 20
dse <- 5
ru <- 3
su <- 5
##################### Amostrador de Gibbs ###
for (i in 1:iter)
{
########## condicional para theta ####
96
V <- rbind(cbind(se
*
diag(3) , matrix(0,3,35)),
cbind(matrix(0,35,3), vu
*
A))
S <- delta2
*
ginv(V)
theta <- ginv(t(W)%
*
%R%
*
%W + S)%
*
%t(W)%
*
%R%
*
%L
var <- delta2
*
ginv(t(W)%
*
%R%
*
%W+ S)
thetaest <- rmvnorm(1,mean = theta,sigma = var, method = "chol")
theta <- t(thetaest)
beta <- theta[1:3]
u <- theta[4:38]
### condicional para variancia de blocos ###
c1 <- (35+2
*
ru)/2
c2 <- (t(u)%
*
%ginv(A)%
*
%u + 2
*
su)/2
vu <- rgamma(1,c1,c2)
vu <- as.real(1/vu)
## par. de aceleramento da convergência (reparametrização) ##
delta2 <- 1/rgamma(1,(n+38+5+2
*
dre)/2,(t(L-W%
*
%theta)%
*
%(L-W%
*
%theta) +
2
*
dse +t(theta)%
*
%ginv(V)%
*
%theta)/2)
delta <- sqrt(delta2)
############## Amostra lambda (GS elemento a elemento) #####
lambda <- rgamma(n,(nu + 1)/2,(nu+(L-W%
*
%theta)^2)/2)
R <- diag(lambda)
################ Amostra nu (MH) ##########################
nuc <- 2
while(nuc < 3)
{
nuc <- rpois(1,nu)
}
pnu <- n
*
((nu/2)
*
log(nu/2)-lgamma(nu/2))+sum(((nu/2)-1)
*
log(lambda)+
(-(nu/2)
*
lambda)) - 2
*
log(1+nu)
pnuc <- n
*
((nuc/2)
*
log(nuc/2)-lgamma(nuc/2))+sum(((nuc/2)-1)
*
log(lambda)+
(-(nuc/2)
*
lambda)) - 2
*
log(1+nuc)
pnu <- exp(pnu)
pnuc <- exp(pnuc)
dnu <- dpois(nu,nu)
dnuc <- dpois(nuc,nu)
num <- pnuc
*
dnu
denom <- pnu
*
dnuc
ifelse((num > denom),(alfa <- 1),(alfa <- as.real(num/denom)))
teste <- runif(1)
if(teste < alfa)
{
nu <- nuc
}
############# Amostragem dos par. de limiar #############
97
nkp <- 0.8
*
nk[2:4]
pnovo <- c(rdiric(1, shape = nkp ))
for(j in 2:4)
{
taunovo[j]<-sum(pnovo[1:(j-1)])
}
lfp <- nkp
*
log(p)
lfpnovo <- nkp
*
log(pnovo)
m <- W%
*
%theta
for(k in 1:n)
{
dp <- rep(1,n)
dp <- sqrt(delta2/lambda)
if(y[k]==1)
{
ftau[k] <- 1
ftaunovo[k]<- 1
}
if(y[k]==2)
{
ftau[k] <- pnorm((tau[2] - m[k])/dp[k]) -
pnorm((tau[1] - m[k])/dp[k])
ftaunovo[k]<- pnorm((taunovo[2] - m[k])/dp[k]) -
pnorm((taunovo[1] - m[k])/dp[k])
}
if(y[k]==3)
{
ftau[k] <- pnorm((tau[3] - m[k])/dp[k]) -
pnorm((tau[2] - m[k])/dp[k])
ftaunovo[k]<- pnorm((taunovo[3] - m[k])/dp[k]) -
pnorm((taunovo[2] - m[k])/dp[k])
}
if(y[k]==4)
{
ftau[k] <- pnorm((tau[4] - m[k])/dp[k]) -
pnorm((tau[3] - m[k])/dp[k])
ftaunovo[k]<- pnorm((taunovo[4] - m[k])/dp[k]) -
pnorm((taunovo[3] - m[k])/dp[k])
}
if(y[k]==5)
{
ftau[k] <- 1
ftaunovo[k]<- 1
}
}
lftau <- log(ftau)
lftaunovo<- log(ftaunovo)
alfa <- exp(sum(lftaunovo-lftau)-sum(lfpnovo-lfp))
teste <- runif(1)
D <- min(1,alfa)
if(teste < D)
{
98
p <- pnovo
tau <- taunovo
################## condicional para L ###################
m <- W%
*
%theta
for(k in 1:n)
{
if(y[k]==1)
{
L[k] <-qnorm(pnorm(tau[1],m[k],sd=dp[k])
*
runif(1),mean=m[k],sd=dp[k])
vero[k]<-pnorm((tau[1] - m[k])/dp[k])
}
if(y[k]==2)
{
L[k] <-qnorm(pnorm(tau[1],m[k],sd=dp[k])+(pnorm(tau[2],m[k],sd=dp[k])
-pnorm(tau[1],m[k],sd=dp[k]))
*
runif(1),mean=m[k],sd=dp[k])
vero[k]<-pnorm((tau[2] - m[k])/dp[k]) - pnorm((tau[1] - m[k])/dp[k])
}
if(y[k]==3)
{
L[k] <-qnorm(pnorm(tau[2],m[k],sd=dp[k])+(pnorm(tau[3],m[k],sd=dp[k])
-pnorm(tau[2],m[k],sd=dp[k]))
*
runif(1),mean=m[k],sd=dp[k])
vero[k]<-pnorm((tau[3] - m[k])/dp[k]) - pnorm((tau[2] - m[k])/dp[k])
}
if(y[k]==4)
{
L[k] <-qnorm(pnorm(tau[3],m[k],sd=dp[k])+(pnorm(tau[4],m[k],sd=dp[k])
-pnorm(tau[3],m[k],sd=dp[k]))
*
runif(1),mean=m[k],sd=dp[k])
vero[k] <-pnorm((tau[4] - m[k])/dp[k]) - pnorm((tau[3] - m[k])/dp[k])
}
if(y[k]==5)
{
L[k] <-qnorm(pnorm(tau[4],m[k],sd=dp[k]) +
(1-pnorm(tau[4],m[k],sd=dp[k]))
*
runif(1),mean=m[k],sd=dp[k])
vero[k]<-(1 - pnorm( - m[k]/dp[k]))
}
}
}
############# Cálculo da verossimilhança #####################
verossim <- prod(vero)
############ Gera a saída no arquivo cadeia.txt ##############
saida <- cbind(t(theta),t(tau),delta2,nu,vu,verossim)
colunas <- length(saida)
if(i > burn && (i-burn)%%jump == 0)
{
write(saida, file="cadeia.txt", ncol=colunas, append=TRUE)
}
}
99
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