5.9 Intervalo/Região de Credibilidade

Seção 3.3 de (Paulino, Turkman, and Murteira 2003)
Seção 8.4 de (S. J. Press 2003)
(Hyndman 1996) e (Hyndman, Einbeck, and Wand 2021)
(Meredith and Kruschke 2022)
(Kumar et al. 2019) e (Martin, Kumar, and Lao 2021) (Obrigado, Pedro Roussos e Cristiano Ferrazzo!)
https://web.stanford.edu/class/archive/cs/cs109/cs109.1218/files/student_drive/8.2.pdf (Obrigado, Guilherme Scherer!)

Exemplo 5.15 Considerando novamente os dados do Exemplo 5.6 e uma priori \(Beta(1,1)\), pode-se calcular as estimativas por intervalo de credibilidade da posteriori \(Beta(10,4)\).

# \theta ~ Beta(1,1), \theta|x ~ Beta(1+s,1+f)
s <- 9
f <- 3

# Intervalo de Credibilidade via quantis (maior amplitude)
qbeta(.025, s+1, f+1)

## [1] 0.4618685

qbeta(.975, s+1, f+1)

## [1] 0.9090796

qbeta(.975, s+1, f+1) - qbeta(.025, s+1, f+1) # Amplitude

## [1] 0.4472111

# Intervalo de Credibilidade via HDI de Meredith and Kruschke (2022)
tst <- rbeta(1e5, s+1, f+1)
(hdint <- HDInterval::hdi(tst))

##     lower     upper 
## 0.4849837 0.9247811 
## attr(,"credMass")
## [1] 0.95

diff(hdint) # Amplitude

##     upper 
## 0.4397974

# Intervalo de Credibilidade via HDI de Hyndman (1996)
hdrcde::hdr.den(tst, prob = 95)

## $hdr
##         [,1]     [,2]
## 95% 0.491742 0.932006
## 
## $mode
## [1] 0.7425832
## 
## $falpha
##        5% 
## 0.7128508

hdrcde::hdr.boxplot(tst)

Exemplo 5.16 Em Python.

import numpy as np
from scipy.stats import beta
import matplotlib.pyplot as plt
from hdintervals import hdi
import hdrcde


# Parâmetros da distribuição Beta
s = 9  # Número de sucessos
f = 3  # Número de fracassos

# Intervalo de Credibilidade via quantis (maior amplitude)
q_lower = beta.ppf(0.025, a=s+1, b=f+1)
q_upper = beta.ppf(0.975, a=s+1, b=f+1)
print(f"Intervalo via quantis: [{q_lower:.4f}, {q_upper:.4f}]")
print(f"Amplitude: {q_upper - q_lower:.4f}")

# Intervalo de Credibilidade via HDI de Meredith and Kruschke (2022)
tst = beta.rvs(a=s+1, b=f+1, size=100000)
hdi_interval = hdi(tst, 0.95)
print(f"\nIntervalo HDI (Meredith and Kruschke): [{hdi_interval[0]:.4f}, {hdi_interval[1]:.4f}]")
print(f"Amplitude: {hdi_interval[1] - hdi_interval[0]:.4f}")

# Intervalo de Credibilidade via HDI de Hyndman (1996)
hdr_interval = hdrcde.hdr(tst, prob=95)
print(f"\nIntervalo HDR (Hyndman): [{hdr_interval[0]:.4f}, {hdr_interval[1]:.4f}]")

# Plotando o boxplot com o intervalo HDR
hdrcde.hdr.boxplot(tst, prob=95)
plt.title('Boxplot com Intervalo HDR')
plt.show()

Exemplo 5.17 A biblioteca ArviZ de (Kumar et al. 2019), desenvolvida para Python, fornece ferramentas para análise exploratória de modelos bayesianos, dentre as quais a função arviz.hdi.

Exercício 5.8 Considere os dados do Exemplo 5.6.
a. Obtenha a posteriori considerando uma priori vaga \(\mathcal{U}(0,1)\) e as prioris (conjugadas) \(Beta(0,0)\), \(Beta(1,1)\) e \(Beta(2,2)\). O que você percebe?
b. Considerando sua priori preferência, obtenha as estimativas pontuais de \(\theta\) sob o prisma bayesiano considerando a média, a mediana e a moda.
c. Considerando sua priori preferência, obtenha o intervalo de credibilidade a partir dos quantis \(P_{0.025}\) e \(P_{0.975}\) e do método da máxima densidade. O que você percebe?

Exemplo 5.18 Em Python.

import numpy as np
from scipy.stats import beta
import matplotlib.pyplot as plt
from hdintervals import hdi

# a) Plotando as curvas das distribuições Beta
x = np.linspace(0, 1, 100)

# Beta(10, 4) - vaga
plt.plot(x, beta.pdf(x, 10, 4), label='Beta(10, 4) - vaga')

# Beta(9, 3) - conjugada com Beta(0, 0)
plt.plot(x, beta.pdf(x, 9, 3), label='Beta(9, 3)')

# Beta(10, 4) - conjugada com Beta(1, 1)
plt.plot(x, beta.pdf(x, 10, 4), label='Beta(10, 4)')

# Beta(11, 5) - conjugada com Beta(2, 2)
plt.plot(x, beta.pdf(x, 11, 5), label='Beta(11, 5)')

plt.xlabel(r'$\theta$')
plt.ylabel('Densidade')
plt.title('Distribuições Beta')
plt.legend()
plt.show()

# b) Calculando a média, mediana e moda da posteriori com priori Beta(2, 2)
alfa0 = 2
beta0 = 2
x = 9
n = 12
alfa1 = alfa0 + x
beta1 = beta0 + n - x

media = alfa1 / (alfa1 + beta1)
mediana = (alfa1 - 1/3) / (alfa1 + beta1 - 2/3)
moda = (alfa1 - 1) / (alfa1 + beta1 - 2)

print(f"Média: {media:.4f}")      # Output: 0.6923
print(f"Mediana: {mediana:.4f}")  # Output: 0.7083
print(f"Moda: {moda:.4f}")        # Output: 0.7500

# c) Calculando o intervalo de credibilidade de 95%
q_lower = beta.ppf(0.025, a=alfa1, b=beta1)
q_upper = beta.ppf(0.975, a=alfa1, b=beta1)
print(f"\nIntervalo via quantis: [{q_lower:.4f}, {q_upper:.4f}]")

# Calculando o HDI
hdi_interval = hdi(beta.rvs(a=alfa1, b=beta1, size=100000), 0.95)
print(f"Intervalo HDI: [{hdi_interval[0]:.4f}, {hdi_interval[1]:.4f}]")

Referências

———. 1996. “Computing and Graphing Highest Density Regions.” The American Statistician 50 (2): 120–26. https://www.tandfonline.com/doi/abs/10.1080/00031305.1996.10474359.

Hyndman, Rob J, Jochen Einbeck, and Matthew P Wand. 2021. hdrcde: Highest Density Regions and Conditional Density Estimation. https://pkg.robjhyndman.com/hdrcde/.

Kumar, Ravin, Colin Carroll, Ari Hartikainen, and Osvaldo Martin. 2019. “ArviZ a Unified Library for Exploratory Analysis of Bayesian Models in Python.” Journal of Open Source Software 4 (33): 1143. https://doi.org/10.21105/joss.01143.

Martin, Osvaldo A, Ravin Kumar, and Junpeng Lao. 2021. Bayesian Modeling and Computation in Python. CRC Press. https://bayesiancomputationbook.com/welcome.html.

Meredith, Mike, and John Kruschke. 2022. HDInterval: Highest (Posterior) Density Intervals. https://CRAN.R-project.org/package=HDInterval.

Paulino, Carlos Daniel Mimoso, Maria Antónia Amaral Turkman, and Bento Murteira. 2003. Estatı́stica Bayesiana. Fundação Calouste Gulbenkian, Lisboa. http://primo-pmtna01.hosted.exlibrisgroup.com/PUC01:PUC01:puc01000334509.

Press, S James. 2003. Subjective and objective Bayesian statistics: Principles, models, and applications, 2nd. edition. John Wiley & Sons. http://primo-pmtna01.hosted.exlibrisgroup.com/PUC01:PUC01:oclc(OCoLC)587388980.