Beta-binomial
lsj
May 31, 2021
O modelo beta-binomial
a distribuição binomial
A distribuição binomial descreve problemas onde os dados são do tipo 0/1, TRUE/FALSE, Detectado/não-Detectado, Cara/Coroa; em geral sucesso/não-sucesso
A distribuição binomial tem um único parâmetro, \(\theta\): a probabilidade de sucesso em cada tentativa
Probabilidade de \(n\) sucessos em \(N\) tentativas independentes: \[ P(n|\theta,N)=\binom Nn \theta^n (1-\theta)^{N-n} = {N! \over n!(N-n)!} \theta^n (1-\theta)^{N-n}\] média e variância: \[ \bar n = \sum_{n=0}^N n P(n) = N\theta ~~~~~~ \sigma^2 = \sum_{n=0}^N (n-\bar n)^2 P(n) = N\theta(1-\theta) \]
Exemplo:
par(mfrow=c(2,2)) # panel de 2 x 2 figuras
# distribuição binomial com p=0.1
N=15 #número de tentativas
n = seq(0,20,1) # número de sucessos - numeros inteiros positivos
y = dbinom(n,N,prob=0.10)
N=30
plot(n,y,main="distribuição binomial (N,p)=(15,0.1)",type="h")
y = dbinom(n,N,prob=0.10)
plot(n,y,main="distribuição binomial (N,p)=(30,0.1)",type="h")
# distribuição binomial com p=0.5
N=15
y = dbinom(n,N,prob=0.5)
plot(n,y,main="distribuição binomial (N,p)=(15,0.5)",type="h")
N=30
y = dbinom(n,N,prob=0.5)
plot(n,y,main="distribuição binomial (N,p)=(30,0.5)",type="h")
Dado um conjunto de dados, vamos fazer um modelo bayesiano para estimar \(\theta\) (\(0 \le \theta \le 1\))
Vamos escrever a verossimilhança como a distribuição binomial \[ P(\theta|n,N) = \binom{N}{n} \theta^n(1-\theta)^{N-n} \]
e o prior como uma função Beta, já que a variável \(\theta\) é definida entre 0 e 1: \[ Beta(\theta|a,b) = \Bigg[ \frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)}\Bigg] \theta^{a-1} (1-\theta)^{b-1} \] \[(a,b ~ > 0) \]
Exemplos da função beta:
# distribuição beta em R: dbeta(x,a,b)
# vamos gerar uma sequência de números entre 0 e 1 separados por 0.01
# vamos amostrar a função beta nesses valores
x=seq(0,1,0.01)
# vamos fazer uma figura com 3 x 2 subfiguras
par(mfrow=c(3,2))
# vamos considerar valores diferentes de a e b e ver a cara da função:
a=10
b=30
y = dbeta(x,a,b)
plot(x,y,type='l',col='blue', main=sprintf(paste('a,b = ',a,b,sep=' ')))
a=1
b=3
y = dbeta(x,a,b)
plot(x,y,type='l',col='blue', main=sprintf(paste('a,b = ',a,b,sep=' ')))
a=1
b=1
y = dbeta(x,a,b)
plot(x,y,type='l',col='blue', main=sprintf(paste('a,b = ',a,b,sep=' ')))
a=3
b=1
y = dbeta(x,a,b)
plot(x,y,type='l',col='blue', main=sprintf(paste('a,b = ',a,b,sep=' ')))
a=30
b=10
y = dbeta(x,a,b)
plot(x,y,type='l',col='blue', main=sprintf(paste('a,b = ',a,b,sep=' ')))
a=10
b=10
y = dbeta(x,a,b)
plot(x,y,type='l',col='blue', main=sprintf(paste('a,b = ',a,b,sep=' ')))
o posterior do modelo beta-binomial
posterior: \[ P(\theta|D) \propto \theta^{n+a-1}(1-\theta)^{N-n+b-1} \] \[ \propto Beta(\theta|n+a,N-n+b) \]
note que, nesse caso, tanto o prior quanto o posterior são funções Beta: Beta é um prior conjugado da distribuição binomial
Exemplos
Exemplo 1: baseado em https://www.r-exercises.com/2017/12/03/basic-bayesian-inference-for-mcmc-techniques-exercises-part-1/
Um modelo de formação de galáxias prevê que em aglomerados de galáxias ricos a maior parte das galáxias luminosas sejam elípticas ou lenticulares (em baixos redshifts). Queremos analisar se esse é realmente o caso em uma amostra de aglomerados. Este problema é uma adaptação do descrito no link acima.
Seja \(\theta\) a fração de galáxias early type (ET: E e S0). O estudo de uma grande simulação revela que em uma amostra de \(N=1000\) galáxias luminosas em aglomerados ricos, \(n=735\) são desse tipo. Vamos classificar as galáxias como \(y=1\) se forem ET e \(y=0\) se não o forem.
Vamos ver, inicialmente, qual é a distribuição de probabilidade dos dados, plotando a verossimilhança:
plot(dbinom(x = seq(1, 1000, 1), size = 1000, prob = 735/1000),
type = "l", lwd = 3, main = "verossimilhança", ylab = "densidade de probabilidades", xlab = "número de galáxias early type", col='blue')
abline(v = 735, col = "red",lwd=2)
Vamos considerar inicialmente um prior uniforme para \(\theta\); ele corresponde a \(Beta(\theta|a=1,b=1)\)
a=1
b=1
y = dbeta(x,a,b)
plot(x,y,type='l',col='blue', main='prior', ylab = "densidade de probabilidades",xlab=expression(theta),lwd=3)
e o posterior será \[ P(\theta|D,a,b) \propto Beta(\theta|a=736,b=266):\]
a=1
b=1
n=735
N=1000
y = dbeta(x,n+a,N-n+b)
plot(x,y,type='l',col='green', main='posterior', ylab = "densidade de probabilidades",xlab=expression(theta),lwd=3)
# moda
# https://stackoverflow.com/questions/2547402/how-to-find-the-statistical-mode
Mode <- function(x) {
ux <- unique(x)
ux[which.max(tabulate(match(x, ux)))]
}
M=100000
post_sample = rbeta(M, shape1 = n+a, shape2 = N-n+b)
Mode(post_sample)## [1] 0.7391225# intervalo de confian,ca de 95%:
quantile(post_sample, probs = c(0.025, 0.975))## 2.5% 97.5%
## 0.7066068 0.7613735Uma outro estudo indica que apenas 60+/-20% das galáxias luminosas são ET. Queremos usar esta informação como prior. Como podemos fazer isso?
A distribuição beta tem dois hiperparâmetros (“shape parameters”) que podem ser obtidos da média e variância da distribuição:
estBetaParams <- function(mean, variance) {
alpha <- ((1 - mean) / variance - 1 / mean) * mean^2
beta <- alpha * (1 / mean - 1)
return(params = list(a = alpha, b = beta))
}Assim, em nosso exemplo o prior de \(\theta\) tem parâmetros:
prior.params <- estBetaParams(mean = 0.6, variance = 0.2^2)
prior.params## $a
## [1] 3
##
## $b
## [1] 2Vamos representar graficamente este prior:
x_grid <- seq(from = 0,to = 1, by = 0.001)
# note que seq deve estar no intervalo de Beta, [0,1]
plot(x = x_grid, dbeta(x = x_grid, shape1 = prior.params$a, shape2 = prior.params$b), lwd = 3, col = 'blue', ylab = "densidade", xlab = expression(theta), type = "l", main = "prior")
E o posterior será:
N <- 1000
n <- 735
# hiper-parâmetros do posterior
alpha_post <- prior.params$a + n
beta_post <- prior.params$b + N - n
#plot do posterior
plot(x = x_grid, dbeta(x = x_grid, shape1 = alpha_post, shape2 = beta_post), lwd = 3, type = "l", ylab = "densidade", xlab = expression(theta), col = "red", main = "posterior e prior")
# prior
lines(x = x_grid, dbeta(x = x_grid, shape1 = prior.params$a, shape2 = prior.params$b), lwd = 3)
legend("topright",c("Posterior","Prior"),col=c("red", "black"),lwd=c(3,3))
# valor de teta da verossimilhança
abline(v = 0.735, col = "green")
# moda
M=100000
post_sample = rbeta(M, shape1 = n+a, shape2 = N-n+b)
Mode(post_sample)## [1] 0.7269337# intervalo de confian,ca de 95%:
post_sample = rbeta(M, shape1 = prior.params$a, shape2 = prior.params$b)
quantile(post_sample, probs = c(0.025, 0.975))## 2.5% 97.5%
## 0.1960180 0.9329706Pode-se notar que, nesse exemplo, a forma do prior afeta muito pouco o posterior. Compare a moda e o intervalo de confiança de \(\theta\) nos dois casos. Se o número de dados fosse bem menor as diferenças seriam maiores!
Exemplo 2: este exemplo é do blog de Rasmus Bååth (ver http://www.sumsar.net/) e mostra um exemplo do problema binomial, ilustrando a evolução do posterior de acordo com o número de observações.
https://www.r-bloggers.com/2018/12/yet-another-visualization-of-the-bayesian-beta-binomial-model/
# This function takes a number of successes and failuers coded as a TRUE/FALSE
# or 0/1 vector. This should be given as the data argument.
# The result is a visualization of the how a Beta-Binomial
# model gradualy learns the underlying proportion of successes
# using this data. The function also returns a sample from the
# posterior distribution that can be further manipulated and inspected.
# The default prior is a Beta(1,1) distribution, but this can be set using the
# prior_prop argument.
# Make sure the packages tidyverse and ggridges are installed, otherwise run:
# install.packages(c("tidyverse", "ggridges"))
# Example usage:
# data <- c(TRUE, FALSE, TRUE, TRUE, FALSE, TRUE, TRUE)
# prop_model(data)
prop_model <- function(data = c(), prior_prop = c(1, 1), n_draws = 10000) {
library(tidyverse)
data <- as.logical(data)
# data_indices decides what densities to plot between the prior and the posterior
# For 20 datapoints and less we're plotting all of them.
data_indices <- round(seq(0, length(data), length.out = min(length(data) + 1, 20)))
# dens_curves will be a data frame with the x & y coordinates for the
# denities to plot where x = proportion_success and y = probability
proportion_success <- c(0, seq(0, 1, length.out = 100), 1)
dens_curves <- map_dfr(data_indices, function(i) {
value <- ifelse(i == 0, "Prior", ifelse(data[i], "Success", "Failure"))
label <- paste0("n=", i)
probability <- dbeta(proportion_success,
prior_prop[1] + sum(data[seq_len(i)]),
prior_prop[2] + sum(!data[seq_len(i)]))
probability <- probability / max(probability)
data_frame(value, label, proportion_success, probability)
})
# Turning label and value into factors with the right ordering for the plot
dens_curves$label <- fct_rev(factor(dens_curves$label, levels = paste0("n=", data_indices )))
dens_curves$value <- factor(dens_curves$value, levels = c("Prior", "Success", "Failure"))
p <- ggplot(dens_curves, aes(x = proportion_success, y = label,
height = probability, fill = value)) +
ggridges::geom_density_ridges(stat="identity", color = "white", alpha = 0.8,
panel_scaling = TRUE, size = 1) +
scale_y_discrete("", expand = c(0.01, 0)) +
scale_x_continuous("Underlying proportion of success") +
scale_fill_manual(values = hcl(120 * 2:0 + 15, 100, 65), name = "", drop = FALSE,
labels = c("Prior ", "Success ", "Failure ")) +
ggtitle(paste0(
"Binomial model - Data: ", sum(data), " successes, " , sum(!data), " failures")) +
theme_light() +
theme(legend.position = "top")
print(p)
# Returning a sample from the posterior distribution that can be further
# manipulated and inspected
posterior_sample <- rbeta(n_draws, prior_prop[1] + sum(data), prior_prop[2] + sum(!data))
invisible(posterior_sample)
}
set.seed(123)
big_data <- sample(c(TRUE, FALSE), prob = c(0.75, 0.25),
size = 100, replace = TRUE)
posterior <- prop_model(big_data)## ── Attaching packages ────────────────────── tidyverse 1.3.0 ──## ✓ ggplot2 3.2.1 ✓ purrr 0.3.3
## ✓ tibble 2.1.3 ✓ dplyr 0.8.3
## ✓ tidyr 1.0.2 ✓ stringr 1.4.0
## ✓ readr 1.3.1 ✓ forcats 0.4.0## ── Conflicts ───────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()## Warning: `data_frame()` is deprecated, use `tibble()`.
## This warning is displayed once per session.
quantile(posterior, c(0.025, 0.5, 0.975))## 2.5% 50% 97.5%
## 0.6360817 0.7272447 0.8082072