A função de verossimilhança liga nossos dados a modelos probabilísticos. Enquanto distribuições de probabilidades atribuem probabilidades a cada observação possível, a função de verossimilhança expressa a plausibilidade de diferentes modelos, dada uma mesma observação ou conjunto de dados.
Os tutoriais a seguir mostram como construir funções de verossimilhança em linguagem computacional, e também como encontrar o valor máximo destas funções com otimização computacional. Vamos fazer também algumas simulações para demonstrar propriedades das funções e dos estimadores de máxima verossimilhança, que são o resultado dessas otimizações. Como sempre, o objetivo é usar recursos computacionais para tornar mais concretos conceitos matemáticos e estatísticos.
Dada uma variável aletória $X$, cujo comportamento pode ser explicado por duas hipóteses: $H_A$ e $H_B$.
A observação $X=x$ é uma evidência em favor de $H_A$ vis-a-vis (face-a-face) $H_B$ se, e somente se,
| $ p_A(x) > p_B(x)$. |
A força de evidência em favor de $H_A$ vis-a-vis $H_B$ é dada pela razão de verossimilhança:
| ${p_A(x)} / {p_b(x)}$ |
Numa ninhada de 5 filhotes de cachorros-do-mato foi observado apenas um macho.
A hipótese H1 estabelece que a probabilidade de nascimento de um filhote do sexo masculino é p=1/3, enquanto que a hipótese H2 estabelece uma probabilidade p=1/2 para filhote macho.
A probabilidade de observar 1 macho em 5 é:
> dbinom(1, 5, 1/3) [1] 0.3292181
> dbinom(1, 5, 1/2) [1] 0.15625
Questão: Qual das duas hipóteses é mais plausível, isto é, mais verossímil? R: A hipótese H1, pois atribui mais probabilidade à observação de um macho em uma ninhada de 5 do que a hipótese H2 ($ p_H1(x) > p_H2(x)$
Questão: Qual a força de evidência em favor da hipótese mais plausível (mais verossímil)?
> dbinom(1, 5, 1/3) / dbinom(1, 5, 1/2) [1] 2.106996
R: Usando a razão de verossimilhanças como medida de força de evidência dos dados, a hipótese H1 é duas vezes mais plausível que a hipótese H2. Ou seja, a observação é uma evidência mais forte a favor da hipótese H1 do que da hipótese H2.
A variável $X$ é uma variável aleatória com função de densidade $f_X(x; \theta)$, onde:
Para que a função de densidade $f_X(\cdot)$ seja utilizada para calcular a probabilidade da observação $X=x$ é necessário que o valor dos parâmetros $\theta$ sejam conhecidos.
Na pesquisa científica o valor dos parâmetros nunca são conhecidos. O que a pesquisa gera são as observações a partir das quais se deseja estimar os parâmetros do modelo.
Assim, a função de densidade pode ser utilizada de uma forma diferente: cria-se com ela uma função dos parâmetros, mantidos os dados constantes.
Essa é a função de verossimilhança, que é a função de densidade com parâmetros desconhecidos e condicionada aos dados, que são conhecidos:
$$\mathcal{L} (\theta | X=x) = f_X(\theta | X = x)$$
Assim, na função de verossimilhança o elemento que varia são os parâmetros ($\theta$), enquanto que o elemento constante é a observação ($X=x$).
Considere que numa ninhada de cachorro-do-mato nasceu 1 macho em 5 filhotes.
Assumindo que o número de machos numa ninhada segue a distribuição binomial, o que essa observação nos diz sobre a probabilidade de nascimento de filhotes machos ($p$)?
Vejamos um gráfico da função de verossimilhança para os diferentes valores de $p$ que essa observação gera:
curve(dbinom(1, 5, p=x), from=0, to=1)
Portanto: O gráfico descreve uma curva de plausibilidade para os valores de probabilidade de nascer macho, dada a observação de 1 macho em 5 filhotes 1).
A função nos apresenta os valores de verossimilhança, sendo razoável assumir que eles medem a plausibilidade de cada valor do parâmetro $p$ numa escala relativa, isto é, em comparação com os demais valores do gráfico.
PONTO de MÁXIMO da FUNÇÃO de VEROSSIMILHANÇA
O valor de máxima verossimilhança (= máxima plausibilidade) pode ser encontrado e colocado no gráfico:
theta = seq(0,1,length=100) y = dbinom(1, 5, p=theta) theta.mle = theta[ y == max(y) ] abline( v = theta.mle, lty=2, col="orange" )
VEROSSIMILHANÇA RELATIVA
Podemos, então, re-escalonar o eixo-y (ordenadas) do gráfico em relação a esse valor máximo:
yrel = y / max(y) plot(theta, yrel, type="l", xlab="Probabilidade p", col="blue") abline( v = theta.mle, lty=2, col="orange" )
Portanto: podemos agora chamar esse gráfico de verossimilhança relativa (plausibilidade relativa) para os valores de probabilidade de nascimento de filhote macho, dada a observação de 1 macho em 5 filhotes. [Observe os valores no eixo-y]
INTERVALO DE VEROSSIMILHANÇA
Imagine que queremos encontrar nesse gráfico um intervalo para todos os valores cujo valor de verossimilhança relativa é de pelo menos 1/8 em relação ao máximo. Isto é, o valor máximo tem verossimilhança relativa de até 8 vezes os demais valores do intervalo.
theta.interv = theta[ max(y)/y <= 8 ] lines(theta.interv, rep(1/8, length(theta.interv)), col="red", lwd=3)
Geralmente, os dados científicos são obtidos na forma de uma amostra da população, sendo que essa amostra é composta de várias observações independentes.
Se foram tomadas 4 observações independentes de uma variável aleatória: $X$ ($\mathbf{X}_4 = \{x_1, x_2, x_3, x_4\}$) e cada uma delas possui uma função de verossimilhança para o parâmetro ($\theta$) do modelo, então a função de verossimilhança da amostra é o produto das funções de verossimilhança das observações independentes:
$$ \mathcal{L} ( \theta | X = \mathbf{X}_4 ) = \mathcal{L} ( \theta | X = x_1 ) \cdot \mathcal{L} ( \theta | X = x_2 ) \cdot \mathcal{L} ( \theta | X = x_3 ) \cdot \mathcal{L} ( \theta | X = x_4 )$$
Generalizando: se uma amostra é composta de $n$ observaçoes independentes, então a verossimilhança da amostra é o produto das verossimilhanças das observações individuais:
$$ \mathcal{L} ( \theta | X = \mathbf{X}_n ) = \mathcal{L} ( \theta | X = x_1 ) \cdot \mathcal{L} ( \theta | X = x_2 ) \cdot \mathcal{L} ( \theta | X = x_3 ) \cdot \ldots \cdot \mathcal{L} ( \theta | X = x_n )$$
$$ \mathcal{L} ( \theta | X = \mathbf{X}_n ) = \prod_{i=1}^n \mathcal{L} ( \theta | X = x_i )$$
Num levantamento em floresta ombrófila densa foram contadas o número de espécies por parcelas de 900m2, obtendo-se os valores: 5, 8, 11 e 9. Como se trata dados de contagem, a distribuição de Poisson é candidata à modelagem desses dados.
Se construírmos uma curva de verossimilhança para cada observação temos:
curve(dpois( 5, lambda=x), 0, 30 , xlab=expression(lambda), ylab="") curve(dpois( 8, lambda=x), 0, 30 , add=TRUE) curve(dpois( 9, lambda=x), 0, 30 , add=TRUE) curve(dpois( 11, lambda=x), 0, 30 , add=TRUE)
Entretanto, nossa amostra é formada pelas quatro observações conjuntamente e precisamos representar uma única curva:
x = seq(0,30, by=0.1) y1 = dpois(5, lambda=x) y2 = dpois(8, lambda=x) y3 = dpois(9, lambda=x) y4 = dpois(11, lambda=x) y = y1 * y2 * y3 * y4 plot(x, y, type="l", col="red")
Em um único passo, teríamos:
curve(dpois(5,x)*dpois(8,x)*dpois(9,x)*dpois(11,x), 0, 30 , xlab=expression(lambda), ylab="Verossimilhança", col="red")
A função de verossimilhança quando aplicada a uma única observação resulta frequentemente (mas não necessariamente) em valores menores do que um. Quando ela é aplicada a uma amostra com muitas observações, o produto das verossimilhanças individuais vai rapidamente em direção ao zero, se tornando um número muito pequeno.
Embora se possa trabalhar com a verossimilhança relativa, a função de verossimilhança é geralmente uma função de difícil tratamento matemático.
Se aplicamos o logaritmo e mudança de sinal à função de verossimilhança relativa, obtemos a função de log-verossimilhança negativa que é uma função de tratamento matemático mais simples e que resulta em valores numéricos de manipulação mais fácil.
Se compararmos os valores nos vetores c(y1,y2,y3,y4) e y, veremos que ao fazermos o produto, obtemos números muito pequenos.
range(c(y1,y2,y3,y4)) #[1] 0.0000000 0.1754674 range(y) #[1] 0.0000000000 0.0001162465
Em grandes amostras, com 100 ou mais observações, esse produto gera números muito pequenos, sendo problemático o seu processamento no computador. Por isso é conveniente utilizarmos a transformação para a função de log-verossimilhança negativa:
log.y = - log(y) log.y2= - (log(y1) + log(y2) + log(y3) + log(y4)) plot(x, log.y, type="l", col="red", xlab=expression(lambda)) lines(x, log.y2, col="blue")
Note que log.y e log.y2 são identicos.
No R, a log-verossimilhança negativa pode ser gerada diretamente das funções de densidade das distribuições utilizando o argumento “log = TRUE” que existe em todas elas. Vejamos uma forma mais prática (e elegante) de se gerar a função de log-verossimilhança negativa da distribuição Poisson para esse conjunto de dados:
library(MASS) # carrega o pacote "MASS" args(dpois) # verifica os argumentos da função "dpois" lpois = function(x,lambda) -dpois(x,lambda,log=TRUE) # cria uma função da log-veros. neg. llikpois = Vectorize( lpois, c("x","lambda") ) # cria uma função "vetorizada" da log-veros. neg. x = c(5, 8, 9, 11) # "x" agora são as observações lambda = seq(0, 30, by=0.1) # "lambda" o vetor com os valores do parâmetro lver.mat = outer( x, lambda, "llikpois" ) # uma matrix a log-veros. neg. de cada observação lver = apply( lver.mat, 2, sum) # um vetor com a log-veros. neg da amostra plot( lambda, lver, type="l", col="red", xlab=expression(lambda) )
Da mesma forma que re-escalonamos a curva de verossimilhança em relação ao seu máximo, podemos re-escalonar essa curva de log-verossimilhança negativa em relação ao seu mínimo:
plot(lambda, lver - min(lver), type="l", col="red", xlab=expression(lambda) )
O ponto que maximiza a verossimilhança é o mesmo que minimiza a log-verossimilhança negativa:
lambda.min = lambda[ min(lver) == lver ] abline(v = lambda.min, lty=2, col="purple")
Note que o valor mais verossímil graficamente para o parâmetro $\lambda$ da Poisson é:
lambda.min #[1] 8.3
que é muito próximo a média amostral, a solução analítica para o estimador de máxima verossimilhança $\lambda$ da Poisson:
mean(c(5,8,9,11)) #[1] 8.25
A diferença que encontramos é apenas por causa da precisão desta solução gráfica. Que é didática mas pouco precisa.
Também podemos estabelecer graficamente um intervalo de verossimilhança para o parâmetro $\lambda$ na curva de log-verossimilhança negativa:
lambda.int = lambda[ lver - min(lver) <= log(8) ] lines(lambda.int, rep(log(8), length(lambda.int)), col="purple", lwd=3 )
O método da Máxima Verossimilhança consiste em realizar de uma maneira mais precisa o que foi feito nos exemplos acima, isto é, encontrar a estimativa dos parâmetros de um modelo que maximiza a função de verossimilhança ou, o que é equivalente, minimizar a função de log-verossimilhança negativa.
As estimativas assim geradas são chamadas de MLE = Maximum Likelihood Estimates.
O R possui duas funções básicas para minimizar expressões matemáticas, ou seja, para buscar os valores dos parâmetros de uma função que resultam no menor valor da função:
A partir destas funções há várias definidas especificamente para minimizar funções de verossimilhança. As duas mais usadas são a
mle do pacote stats4, e a mle2, do pacote bbmle2).
Como o pacote bbmle foi feito por um ecólogo, para acompanhar o ótimo livro sobre modelos que usamos como referência básica, usaremos a mle2 em nossos tutoriais.
A função que deve ser fornecida para a “mle2” é semelhante às funções apresentadas acima, mas nesse caso ela deve ser uma função que soma a log-verossimilhança negativa para cada valor dos parâmetros para todas observações na amostra.
Nesse exemplo utilizaremos os dados de número de plântulas do palmiteiro juçara (Euterpe edulis) encontradas em parcelas de regeneração natural. Vamos ajustar a distribuição Poisson ao dados, o que siginifica buscar a estimativa de máxima verossimilhança do único parâmetro da distribuição Poisson 3).
euplant = c(5, 23, 23, 17, 0, 4, 1, 0, 0, 1, 0, 2, 26, 65, 34, 14, 18, 13, 19, 7) length(euplant) mean(euplant)
É necessário criar a função de log-verossimilhança negativa da amostra para distribuição Poisson:
nvl = function(lambda) { -sum( dpois( euplant, lambda, log=TRUE ) ) }
Vamos encontrar a estimativa de máxima verossimilhança (MLE) do parâmetro $\lambda$:
library(bbmle) # Primeiro carrega o pacote "bbmle" euplant.mle = mle2(nvl, start=list(lambda=10) ) # gera um objeto "mle" summary(euplant.mle) # resumo do objeto "mle" logLik(euplant.mle) # valor da log-verossimilhança negativa do objeto
Verifique se a MLE gerada pela função “mle” de fato corresponde à estimativa esperada teoricamente:
coef(euplant.mle) mean(euplant)
Para visualizar a curva da log-verossimilhança negativa use a função plotprofmle, do pacote sads. Instale e carregue o pacote e então será possível visualizar facilmente a curva:
library(sads) plotprofmle( profile(euplant.mle) )
Os estimadores de máxima verossimilhança têm várias propriedades de um estimador ideal. Para demonstrar estas propriedades dos MLEs é necessário lançar mão de um expediente de simulação.
A simulação consiste nos seguintes passos:
A consistência é a propriedade de convergir em probabilidade para o valor do parâmetro com o aumento da amostra. Para verificar a consistência, utilizaremos o MLE do desvio padrão da Gaussiana, que é sabidamente enviesado:
$$ \widehat{s} = \left[ \frac{1}{n} \sum_{i=1}^n (x_i - \bar{x})^2 \right] ^{(1/2)} $$
Primeiramente, criaremos uma função para obter esta MLE:
sd.mle = function(x) sqrt( (1/length(x)) * sum( (x - mean(x))^2 ) )
Esse estimador será aplicado a amostras da distribuição Gaussiana com média 100 e desvio padrão 5, de tamanhos crescentes.
Gerando a distribuição das MLE para amostras de tamanho 10:
samp10 = sort(rep(1:1000, 10)) # identificador das 1000 amostras de tamanho 10 vals10 = rnorm(1000*10, mean=10, sd=5) # valores das observaçoes dens10 = density( tapply(vals10, samp10, sd.mle) ) # obtendo a MLE para cada amostra e a densidade empírica dens10$y = dens10$y/ max(dens10$y) # convertendo o valor de densidade para valor relativo ao max. plot( dens10, type="l", xlab="Valor Estimado", ylab="Densidade" ,main="") # gráfico da distribuição das MLE abline(v = 5, col="red", lwd=2 ) # posição do valor do parâmetro no gráfico
Note que a distribuição das MLE das amostras de tamanho 10 não é simétrica e não parece centrada no valor do parâmetro.
Gerando a distribuição das MLE para amostras de tamanho 100:
samp100 = sort(rep(1:1000, 100)) vals100 = rnorm(1000*100, mean=100, sd=5) dens100 = density( tapply(vals100, samp100, sd.mle) ) dens100$y = dens100$y/ max(dens100$y) lines( dens100, col="blue" )
Note que as MLE das amostras de tamanho 100 já tem distribuição simétrica e o ponto central se aproxima bastante do valor do parâmetro. Mais precisamente, a probabilidade de que probabilidade das estimativas estarem a uma certa distância do valor real fica cada vez menor com o aumento da amostra.
A eficiência assintótica afirma que os MLE terão a menor variância possível à medida que o tamanho de amostra aumenta.
Vamos comparar a variâncias dos MLEs do desvio padrão obtidos acima (nossa função “sd.mle”) contra os valores do estimador tradicional, calculados para as mesmas amostras (a função sd do R):
## Amostras de tamanho 10 var( tapply(vals10, samp10, sd.mle) ) var( tapply(vals10, samp10, sd) ) ## Amostras de tamanho 100 var( tapply(vals100, samp100, sd.mle) ) var( tapply(vals100, samp100, sd) )
Esta é a propriedade de distribuição assintótica que garante que a distribuição dos MLE 4) aproximam-se de uma distribuição normal com o aumento do tamanho da amostra. Pode ser entendida como uma extensão de Teorema Central do Limite para todos os MLEs.
A melhor forma de investigar a aproximação assintótica da distribuição Gaussiana é através do gráfico Quantil-Quantil:
par(mfrow=c(2,1)) qqnorm( tapply(vals10, samp10, sd.mle), main="MLE do Desvio Padrão, N=10" ) qqline( tapply(vals10, samp10, sd.mle) ) qqnorm( tapply(vals100, samp100, sd.mle), main="MLE do Desvio Padrão, N=100") qqline( tapply(vals100, samp100, sd.mle) ) par(mfrow=c(1,1))
A propriedade de invariância dos MLE é de difícil de comprovação por simulação, sendo mais fácil exemplificar com estimadores que não têm essa propriedade.
O estimador tradicional da variância tem a propriedade de não ser viciado, independentemente do tamanho da amostra.
$$ \widehat{s^2} = \frac{ \sum_{i=1}^n (x_i - \bar{x})^2 }{ n - 1 }$$
Isso pode ser comprovado com uma simulação:
var.trad10 = tapply( vals10, samp10, var) # estimativas tradicionais da variância de amostras de tamanho 10 plot( density(var.trad10), type="l" ) # distribuição das estimativas abline( v = mean(var.trad10) ) # média das estimativas abline( v = 5^2, col="red" ) # valor do parâmetro
No entanto, este não é um estimador de máxima verossimilhança. Por isso a propriedade de não ser viciado não se mantém quando a variância é transformada em desvio padrão , que é uma transformação monotônica:
$$ \widehat{s} = \sqrt{ \widehat{s^2}} \ = \ \sqrt{ \frac{ \sum_{i=1}^n (x_i - \bar{x})^2 }{ n - 1 } }$$
Verifique isto repetindo a avaliação de viés feita antes para a variância:
sd.trad10 = tapply(vals10, samp10, sd ) plot(density(sd.trad10), type="l") abline( v = mean(sd.trad10) ) abline( v = 5, col="red" )
Quando temos mais de parâmetro no modelo estatístico, a função de verossimilhança se torna uma função multivariada e, consequentemente, ela só pode ser representada graficamente por superfícies.
No caso de dois parâmetros, como a distribuição Gaussiana, a função de verossimilhança é bivariada, sendo possível construir gráficos de superfície, mas com três ou mais parâmetros a construção de gráficos e sua interpretação se torna complexa.
dados = rnorm(1000, mean=150, sd=20 ) # gerando os dados m = seq(125, 175, length=50) # vetor de médias s = seq(10, 60, length=50) # vetor de desvio padrão lgauss = function(m, s, dados) -sum(stats::dnorm(dados, mean=m, sd=s, log=T)) # função de log-ver. neg. llikgauss = Vectorize(lgauss, c("m","s")) # vetorização da função sup.mat = outer(m, s, llikgauss, dados) # cálculo da superfície contour(m, s, sup.mat, xlab="Média", ylab="Desvio Padrão", nlevels=20) # gráfico de contorno abline(v=150, col="red", lty=2) # linha da média abline(h=20, col="red", lty=2) # linha do desvio padrão persp(m, s, sup.mat, theta=135, phi=10, col="purple") # gráfico de perspectiva
Copaifera langsdorffii é uma espécie árborea que ocorre em várias formações florestais no Estado de São Paulo, sendo muito abundante no cerradão. O arquivo copaifera.csv apresenta o número de árvores em 64 parcelas na Estação Ecológica de Assis.
copa = read.csv("copaifera.csv",header=T) mean(copa$narv) var(copa$narv)
A distribuição espacial da Copaifera é em geral agregada, de forma que o número de árvores por parcela pode ser descrita com uma distribuição Binomial Negativa com parâmetros:
lnbinom.x = function(k, mu, x) -sum(dnbinom(x, size=k, mu=mu, log=T)) # função de log-veros. neg. lnbinom.vec = Vectorize( lnbinom.x, c("k","mu") ) # função vectorizada k.var = seq(0.1, 10, length=50) # vetor do parâmetro k mu.var = seq(10, 80, length=50) # vetor do parâmetro mu sup.mat = outer( k.var, mu.var, lnbinom.vec, x=copa$narv ) # matriz com a log-veros. neg. contour(k.var, mu.var, sup.mat, nlevels=20, xlab="k", ylab=expression(mu)) # gráfico de contorno persp(k.var, mu.var, sup.mat, theta=-35, phi=20, col="lightblue") # gráfico de perspectiva
Quando o modelo estatístico tem muitos parâmetros, soma-se ao problema de interpretar uma superfície multivariada a questão de alguns dos parâmetros do modelo serem incovenientes, isto é, não possuem uma interpretação biológica clara, mas são necessários ao modelo.
Há duas formas de contornar esse problema que permitem investigar somente o comportamento dos parâmetros de interesse um-a-um, sem a necessidade de avaliar uma superfície:
A verossimilhança perfilhada é uma forma mais apropriada de estudar a superfície de verossimilhança, pois representa de modo mais fidedigno a superfície de verossimilhança na vizinhança das MLE.
Mas em algumas situações, quando não se tem à disposição os dados primários, a verossimilhança estimada pode ser a única alternativa para se analisar os dados.
Vejamos alguns exemplos:
Num inventário em floresta plantada de Pinus foram medidas 52 parcelas, tendo sido obtida a produção média de 243 $m^3\, ha^{-1}$, com desvio padrão de 28 $m^3\, ha^{-1}$. Qual o intervalo de verossimilhança (razão de 8) dessas estimativa da produção 5) ?
Pelo Teorema Central do Limite, sabemos que a média amostral segue distribuição Gaussiana, com média igual à média populacional e desvio padrão igual ao desvio padrão populacional dividido pela raiz quadrada do tamanho da amostra.
Portanto, podemos assumir que a produção média relatada é uma observação de uma distribuição Gaussiana. Como só dispomos dessa informação poderemos avaliar apenas a variação da média, deixando o desvio padrão constante (igual ao relatado):
prod.var = seq(230, 260, length=100) # intervalo de variação da produção média lvn.prod = -dnorm(243, mean=prod.var, sd=28/sqrt(52), log=T) # valores de log-veros. negativa lvn.prod = lvn.prod - min(lvn.prod) # valores de log-veros. negativa RELATIVA plot(prod.var, lvn.prod, type="l") # gráfico log-veros. neg. relativa abline(h=log(8), lty=2, col="red") # intervalo de verossimilhança
Numa floresta planta de eucalipto foram encontradas 7 árvores bifurcadas numa amostra aleatória de 150. Deseja-se a taxa de bifurcação e um intervalo de verossimilhança.
Neste caso a estimativa é o parâmetro $p$ de uma binomial, sendo o outro parâmetro (número de tentativas, $N$) conhecido.
p.var = seq(0.01, 0.10, length=100) lvn.binom = -dbinom(7, 150, prob=p.var, log=T) lvn.binom = lvn.binom - min(lvn.binom) plot( p.var, lvn.binom, type="l", col="red", xlab="Prob. Bifurcação") p.mle = p.var[lvn.binom == min (lvn.binom)] abline(v=p.mle, lty=2, col="blue") abline(h=log(8), lty=2, col="blue")
O arquivo diam.csv apresenta dados do diâmetro à altura do peito (DAP) de árvores de uma parcela em floresta nativa no município de Bom Jardim, MA.
Frequentemente se assume que a distribuição dos DAP das árvores em floresta tropical segue uma distribuição exponencial. Existe, entretanto, a distribuição Weibull, bem mais geral, da qual a distribuição exponencial é um caso particular. Quando o parâmetro da forma (shape) da Weibull é 1, ela se reduz a uma exponencial.
Assim, podemos verificar se a distribuição exponencial é de fato um bom modelo para os dados acima, comparando-a com a Weibull.
Inicialmente é necessário ler os dados do arquivo
dap = read.csv("diam.csv",header=T)$dap
Para ajustar a distribuição Weibull aos dados é necessário considerar que só foram medidas as árvores com DAP a partir de 14.3 cm (45cm de circunferência) e, portanto, o valor 14.3 representa um valor mínimo conhecido.
dap.a = dap - 14.3 # 14.3 é o DAP mínimo de medição lweibull = function(esc, forma){ -sum( dweibull(dap.a, shape=forma, scale=esc, log=T) ) } # função log-ver.neg. da Weibull dap.mle = mle2( lweibull, start=list(esc=10, forma=1) ) # obtendo as MLE summary(dap.mle) # resumo das MLE dap.mle.p = profile( dap.mle ) # perfilhando as MLE par(mfrow=c(1,2)) # preparando para 2 gráficos plotprofmle( dap.mle.p ) # gráficos da veros. perfilhada
Questão: é razoável assumir que os dados podem ser modelados por uma distribuição exponencial?
Para obtermos a curva de verossimilhança perfilhada para o modelo de Binomial Negativa nos dados de número de árvores de Copaifera langsdorffii devemos re-escrever a função de log-verossimilhança negativa e ajustar o modelo com a função mle2.
lnbinom = function(k, mu) -sum(dnbinom(narv, size=k, mu=mu, log=T)) narv = copa$narv copa.mle = mle2( lnbinom, start=list(k=2, mu=32) ) summary(copa.mle) par(mfrow=c(1,2)) plotprofmle( profile(copa.mle) ) par(mfrow=c(1,1)) contour(k.var, mu.var, sup.mat, nlevels=20, xlab="k", ylab=expression(mu)) abline(v=coef(copa.mle)["k"], col="red") abline(h=coef(copa.mle)["mu"], col="red")
Faça os exercícios 203.1 a 203.3 no sistema notaR.
Proponha questões para discutirmos em nosso subfórum. Criamos um subfórum para isso aqui.
bbmle. Ótimo roteiro para aplicação do método de máxima verossimilhança em R.Nossos notebooks interativos em Python demonstrando os passos para a dedução de alguns estimadores de máxima verossimilhança (MLEs).