Discrepância entre binomial e beta em R?

Question

Estou recebendo um resultado que não posso explicar ao usar a distribuição beta.

Obtive um resultado que veio de uma distribuição binomial: 2 sucessos em 6 tentativas. Eu pensaria que o estimador de máxima verossimilhança para p seria 2/6 = 0,33?

dbinom(0:6, 6, 0.33)
[1] 0.090458382 0.267324771 0.329168562 0.216170399 0.079853991 0.015732428 0.001291468

Mas, quando uso a distribuição beta, o ponto mais alto que obtenho é 0,25:

beta_df <- data.frame(PROB = seq(0, 1, 0.01), HEIGHT = dbeta(seq(0, 1, 0.01), 2, 4))
beta_df[which.max(beta_df$HEIGHT),] beta_df[which.max(beta_df$HEIGHT),]
   PROB   HEIGHT
26 0.25 2.109375

Não consigo entender isso ... estou interpretando mal os resultados ou chamando alguma dessas funções incorretamente? Obrigado :)

Answer 1

Por que você esperaria ver resultados semelhantes? Essas são distribuições diferentes, usadas para modelar coisas completamente diferentes. O primeiro é uma distribuição discreta, o segundo é uma distribuição contínua. Respondendo à sua pergunta, você está procurando o modo das distribuições. O modo de distribuição beta é$\frac{\alpha - 1}{\alpha + \beta - 2}$, tão exatamente $0.25$ no caso dos valores que você forneceu.

Em relação ao comentário, no caso binomial, você estava maximizando a probabilidade sozinho. Ao usar o modelo beta-binomial bayesiano , ao maximizá-lo, você está considerando também o modelo anterior

$$ \hat p = \operatorname{arg\,max} \; \underbrace{p(X|\theta)}_\text{likelihood}\,\underbrace{p(\theta)}_\text{prior} $$

portanto, a escolha do prior afetaria o resultado. Ao usar$\alpha=\beta=0$no anterior, este é um prior de Haldane impróprio que tem toda a massa de probabilidade sobre os valores$0$ e $1$(veja a imagem abaixo emprestada deste site ).

Especialmente quando o tamanho da amostra é pequeno , o anterior impactaria o resultado. Nesse caso, ele arrastará a massa de probabilidade para os extremos. Para obter um resultado comparável ao MLE, você pode escolher o uniforme antes com$\alpha=\beta=1$.

Answer 2

Você está usando o parâmetro errado para a distribuição beta. Se você tiver um experimento binomial com$n$ sucessos e $m$ falhas, você deve usar o $beta(n+1,m+1)$distribuição. O motivo é que você está basicamente usando um$beta(1,1)$ (uniforme) anterior, que você deve adicionar à distribuição (se você usar um $beta(a,b)$ antes, em vez disso, você obtém $beta(n+a,m+b)$)

Então, se você tentar isso

beta_df <- data.frame(PROB = seq(0, 1, 0.01), HEIGHT = dbeta(seq(0, 1, 0.01), 3, 5))
beta_df[which.max(beta_df$HEIGHT),]

você obtém o resultado correto, ou seja, $ 0,33 $ .

EDITAR (mais matemática):

Portanto, a probabilidade de $ n $ sucessos e $ m $ falhas é (até uma constante multiplicativa):

$ L (\ theta | \, m, n) \ sim \ theta ^ n (1- \ theta) ^ m $

Mas a densidade beta é dada como

$ f_ {a, b} (\ theta) \ sim \ theta ^ {(a-1)} (1- \ theta) ^ {(b-1)} $ .

Portanto, se você combinar $ a-1 = n $ e $ b-1 = m $, você obterá $ a = n + 1 $ e $ b = m + 1 $ , como deveria.