5.8 Intervalo/Região de Credibilidade

Many statistical methods involve summarizing a probability distribution by a region of the sample space covering a specified probability. (Hyndman 1996, 120)

Exemplo 5.32 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)\) com base em (Meredith and Kruschke 2022), (Hyndman, Einbeck, and Wand 2021) e (Hyndman 1996).

# \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, Meredith and Kruschke (2022)
tst <- rbeta(1e5, s+1, f+1)
(hdint <- HDInterval::hdi(tst))
##     lower     upper 
## 0.4875253 0.9259841 
## attr(,"credMass")
## [1] 0.95
diff(hdint) # Amplitude
##     upper 
## 0.4384587
# Intervalo de Credibilidade via HDI de Hyndman (1996)
hdrcde::hdr.den(tst, prob = 95)

## $hdr
##          [,1]     [,2]
## 95% 0.4896754 0.928572
## 
## $mode
## [1] 0.748938
## 
## $falpha
##        5% 
## 0.6555316
hdrcde::hdr.boxplot(tst)

Exemplo 5.33 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().

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

# Parâmetros da distribuição Beta
s = 9
f = 3

# Intervalo de Credibilidade via quantis (maior amplitude)
lower_quantile = stats.beta.ppf(0.025, s + 1, f + 1)
upper_quantile = stats.beta.ppf(0.975, s + 1, f + 1)
amplitude_quantile = upper_quantile - lower_quantile

print(f"Intervalo de Credibilidade via quantis: ({lower_quantile},{upper_quantile})")
print(f"Amplitude: {amplitude_quantile}")

# Intervalo de Credibilidade via HDI, Meredith and Kruschke (2022)
tst = np.random.beta(s + 1, f + 1, size=100000)
hdint = hdi(tst, hdi_prob=0.95)
amplitude_hdi = hdint[1] - hdint[0]

print(f"Intervalo de Credibilidade via HDI: ({hdint[0]},{hdint[1]})")
print(f"Amplitude: {amplitude_hdi}")

# Intervalo de Credibilidade via HDR de Hyndman (1996)
# Para isso, podemos usar a função `kdeplot` do ArviZ para
# visualizar a densidade e calcular o HDR manualmente ou usar 
# outras bibliotecas como `seaborn` para visualização.

# Exemplo de visualização da densidade com ArviZ
import arviz as az

az.plot_kde(tst, hdi_prob=0.95)
plt.title('Densidade com HDI de 95%')
plt.show()

# Para o HDR boxplot, não há uma função direta equivalente em
# ArviZ, mas podemos usar a função `hdrboxplot` da biblioteca 
# `hdrcde` em R ou implementar manualmente em Python.

Exercício 5.25 Considere os dados do Exemplo 5.6.

  1. 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?
  2. Considerando sua priori preferência, obtenha as estimativas pontuais de \(\theta\) sob o prisma bayesiano considerando a média, a mediana e a moda.
  3. 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.34 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) 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}]")

Exercício 5.26 Veja

  1. Seção 3.3 de (Paulino, Turkman, and Murteira 2003).
  2. Seção 8.4 de (S. J. Press 2003).
  3. (Kumar et al. 2019) e (O. A. Martin, Kumar, and Lao 2021)37.
  4. Seção 8.2 Credible Intervals de (Tsun 2020)38.

References

———. 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. Boca Raton: 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.
Tsun, Alex. 2020. “Probability & Statistics with Applications to Computing.” E-book. https://www.alextsun.com/files/Prob_Stat_for_CS_Book.pdf.

  1. Agradeço a Pedro Roussos e Cristiano Ferrazzo.↩︎

  2. Agradeço a Guilherme Scherer.↩︎