30. Distribuzioni coniugate#
Obiettivo di questo capitolo è fornire un esempio di derivazione della distribuzione a posteriori scegliendo quale distribuzione a priori una distribuzione coniugata. Esamineremo qui il lo schema beta-binomiale.
import pandas as pd
from matplotlib import pyplot as plt
import seaborn as sns
import arviz as az
from scipy import stats
from scipy import integrate
import numpy as np
import pymc as pm
from scipy.stats import beta
%matplotlib inline
%config InlineBackend.figure_format = 'retina'
# Initialize random number generator
RANDOM_SEED = 8927
rng = np.random.default_rng(RANDOM_SEED)
plt.style.use("bmh")
plt.rcParams["figure.figsize"] = [10, 6]
plt.rcParams["figure.dpi"] = 100
plt.rcParams["figure.facecolor"] = "white"
sns.set_theme(palette="colorblind")
%load_ext autoreload
%autoreload 2
%config InlineBackend.figure_format = "svg"
30.1. Derivazione della distribuzione a posteriori#
Le distribuzioni a priori coniugate sono una categoria speciale di distribuzioni probabilistiche che hanno una proprietà importante: se la distribuzione a priori appartiene a questa categoria, allora la distribuzione a posteriori appartiene anche alla stessa categoria, mantenendo la stessa forma funzionale. Ciò significa che l’aggiornamento delle nostre credenze sul parametro di interesse si riduce semplicemente alla modifica dei parametri della distribuzione a priori. Ad esempio, se utilizziamo una distribuzione a priori Beta e la verosimiglianza è binomiale, allora la distribuzione a posteriori sarà anche una distribuzione Beta.
Anche se le distribuzioni a priori coniugate sono la scelta migliore in termini matematici poiché consentono di calcolare la distribuzione a posteriori in modo analitico senza la necessità di calcoli complessi, le tecniche di inferenza bayesiana moderne consentono di utilizzare qualsiasi distribuzione a priori e non solo quelle coniugate. Tuttavia, le distribuzioni a priori coniugate rimangono uno strumento utile nell’insegnamento dell’inferenza bayesiana. In questo capitolo, esploreremo il modello beta-binomiale, che viene utilizzato per l’inferenza sulla proporzione e che si basa sulla distribuzione a priori Beta e sulla verosimiglianza binomiale.
30.2. Lo schema beta-binomiale#
La distribuzione Beta è una distribuzione di probabilità utilizzata per descrivere la variabilità di una variabile casuale che assume valori nell’intervallo [0,1]. La forma della distribuzione è determinata da due parametri, chiamati
Se scegliamo la distribuzione Beta, la distribuzione a priori è
Ometto qui il fattore di normalizzazione di cui non ho bisogno perché la normalizzazione verrà eseguita dopo l’aggiornamento bayesiano.
Per una proporzione, la verosimiglianza è data dalla distribuzione Binomiale:
Nuovamente, ometto il fattore di normalizzazione.
Per calcolare la distribuzione a posteriori, si moltiplica la funzione nucleo a priori Beta per la funzione nucleo della verosimiglianza Binomiale:
La formula risultante rappresenta la distribuzione Beta non normalizzata con i parametri
L’esempio appena descritto rappresenta un caso di analisi coniugata. In particolare, la combinazione della funzione di verosimiglianza Binomiale e della distribuzione a priori Beta è nota come il caso coniugato “beta-binomiale”, ed è regolato dal seguente teorema.
Teorema
Supponiamo di avere una funzione di verosimiglianza
In altre parole, le informazioni a priori rappresentate dalla distribuzione Beta vengono integrate con le informazioni contenute nei dati osservati rappresentati dalla funzione di verosimiglianza Binomiale per produrre una stima più accurata della distribuzione di probabilità a posteriori del parametro
30.2.1. Un esempio concreto#
Per fare un esempio concreto, consideriamo i dati di Zetsche et al. [ZBR19]. In uno studio clinico, 23 partecipanti su 30 hanno riportato aspettative future negative, mentre i restanti 7 hanno riportato aspettative future positive. Indichiamo con
In questo caso, i dati osservati (
30.2.2. La distribuzione a priori#
La distribuzione Beta ci permette di esprimere le nostre iniziali credenze su
def plot_beta(alpha, beta, mean=False, mode=False):
theta = np.linspace(0, 1, 100)
p_theta_given_y = stats.beta.pdf(theta, alpha, beta)
plt.plot(theta, p_theta_given_y, lw=2, color="k")
plt.xlabel(r"$\theta$")
plt.ylabel("$f(\\theta)$")
# Plot the Beta(45, 55) prior
plot_beta(4, 4)
Un sommario della distribuzione
def summarize_beta(alpha, beta):
r"""Summarize a Beta Model for \eqn{\pi}
@param alpha,beta positive shape parameters of the Beta model
Return Pandas Series with summary
"""
mean = alpha / (alpha + beta)
var = alpha * beta / ((alpha + beta) ** 2 * (alpha + beta + 1))
sd = np.sqrt(var)
if alpha < 1 and beta < 1:
mode = "0 and 1"
elif alpha <= 1 and beta > 1:
mode = 0
elif alpha > 1 and beta < 1:
mode = 1
else:
mode = (alpha - 1) / (alpha + beta - 2)
return pd.Series({"mean": mean, "mode": mode, "var": var, "sd": sd})
summarize_beta(4, 4)
mean 0.500000
mode 0.500000
var 0.027778
sd 0.166667
dtype: float64
Possiamo quantificare la nostra incertezza calcolando, con un grado di fiducia del 95%, la regione nella quale, in base a tale credenza a priori, si trova il valore del parametro. Per ottenere tale intervallo di credibilità a priori, usiamo la funzione beta.ppf
di scipy.stats
.
li = beta.ppf(0.025, 4, 4)
ls = beta.ppf(0.975, 4, 4)
list([li, ls])
[0.184051567640083, 0.8159484323599169]
Se poniamo
plot_beta(10, 10)
Tuttavia, in questo caso la nostra certezza a priori sul valore del parametro è maggiore, come indicato dall’intervallo di ordine 0.95.
li = beta.ppf(0.025, 10, 10)
ls = beta.ppf(0.975, 10, 10)
list([li, ls])
[0.2886432479169988, 0.7113567520830011]
Quale distribuzione a priori dobbiamo scegliere? In un problema concreto di analisi dei dati, la scelta della distribuzione a priori dipende dalle credenze a priori che vogliamo includere nell’analisi dei dati. Se non abbiamo alcuna informazione a priori, allora è possibile usare
Nella discussione presente, quale distribuzione a priori useremo una
plot_beta(2, 10)
La distribuzione a priori rappresentata dalla
30.2.3. La distribuzione a posteriori#
Una volta scelta una distribuzione a priori
Essendo
Rappresentiamo qui sotto graficamente l’aggiornamento bayesiano beta-binomiale per i dati di Zetsche et al. [ZBR19] nel caso di una distribuzione a priori
def plot_beta_binomial(
alpha, beta, y=None, n=None, prior=True, likelihood=True, posterior=True
) -> None:
"""Plot a Beta-Binomial Bayesian Model
@param alpha,beta positive shape parameters of the prior Beta model
@param y observed number of successes
@param n observed number of trials
@param prior a logical value indicating whether the prior model should be plotted
@param likelihood a logical value indicating whether the scaled likelihood should be plotted
@param posterior a logical value indicating whether posterior model should be plotted
"""
if y is None or n is None:
print("Warning: to visualize the posterior specify function parameters y and n")
θ = np.linspace(0, 1, 100)
p_theta_given_y = stats.beta.pdf(θ, alpha, beta)
plt.fill_between(θ, p_theta_given_y, lw=4, color="k", label="prior", alpha=0.2)
alpha_post = alpha + y
beta_post = beta + n - y
p_theta_given_y_post = stats.beta.pdf(θ, alpha_post, beta_post)
plt.fill_between(
θ, p_theta_given_y_post, lw=4, color="r", label="posterior", alpha=0.2
)
likelihood = stats.binom.pmf(y, n, θ)
scale_factor = integrate.simpson(likelihood, θ)
plt.fill_between(
θ,
likelihood / scale_factor,
lw=4,
color="b",
label="likelihood scaled",
alpha=0.2,
)
plt.xlabel(r"$\theta$")
plt.ylabel("density")
plt.legend()
plot_beta_binomial(alpha=2, beta=10, y=23, n=30)
30.2.4. La ricerca sull’obbedienza di Milgram#
Consideriamo un altro esempio relativo alla ricerca di Stanley Milgram discussa da Johnson et al. [JOD22]. Nel 1963, Stanley Milgram presentò una ricerca sulla propensione delle persone a obbedire agli ordini di figure di autorità, anche quando tali ordini possono danneggiare altre persone Milgram [Mil63]. Nell’articolo, Milgram descrive lo studio come
consist[ing] of ordering a naive subject to administer electric shock to a victim. A simulated shock generator is used, with 30 clearly marked voltage levels that range from IS to 450 volts. The instrument bears verbal designations that range from Slight Shock to Danger: Severe Shock. The responses of the victim, who is a trained confederate of the experimenter, are standardized. The orders to administer shocks are given to the naive subject in the context of a “learning experiment” ostensibly set up to study the effects of punishment on memory. As the experiment proceeds the naive subject is commanded to administer increasingly more intense shocks to the victim, even to the point of reaching the level marked Danger: Severe Shock.
All’insaputa del partecipante, gli shock elettrici erano falsi e l’attore stava solo fingendo di provare il dolore dello shock.
Johnson et al. [JOD22] fanno inferenza sui risultati dello studio di Milgram mediante il modello Beta-Binomiale. Il parametro di interesse è la probabilità
Sia
Il processo di aggiornamento bayesiano è descritto dalla figura ottenuta con la funzione plot_beta_binomial()
.
plot_beta_binomial(alpha=1, beta=10, y=26, n=40)
30.3. Inferenza bayesiana con distribuzioni a priori continue#
L’inferenza bayesiana sulla proporzione
Il primo tipo (verifica di ipotesi bayesiana) riguarda problemi in cui siamo interessati a valutare la plausibilità che il parametro
Il secondo tipo di inferenza (intervalli di credibilità) riguarda invece la stima dell’intervallo che contiene il parametro
30.3.1. Verifica di ipotesi bayesiana#
Nell’esempio relativo di dati di Milgram, la nostra credenza a posteriori relativa al parametro
Per rispondere a questa domanda, possiamo utilizzare la distribuzione a posteriori Beta(27, 24) e calcolare la probabilità che
1 - stats.beta.cdf(0.5, 27, 24)
0.6640944831173172
30.3.2. Intervalli di credibilità#
Un secondo tipo di inferenza bayesiana è quella che ci porta a costruire gli intervalli di credibilità. Un intervallo di credibilità di ordine
Dato che conosciamo la funzione a posteriori, possiamo semplicemente calcolare i quantili, al livello di probabilità desiderato, per calcolare l’intervallo di credibilità che lascia la stessa probabilità nelle due code.
li = stats.beta.ppf(0.025, 27, 24)
ls = stats.beta.ppf(0.975, 27, 24)
list([li, ls])
[0.3932419761169223, 0.66339490839654]
In alternativa, possiamo calcolare la regione con la più alta densità a posteriori, ovvero la regione (non è necessariamente un intervallo) più corta che contiene la frazione
Questo risultato può essere trovato usando la funzione hdi
del modulo arviz
. Come input della funzione dobbiamo fornire un campione di valori dalla distribuzione a posteriori.
nsim = 1000000
theta_samples = np.random.beta(27, 24, size=nsim)
az.hdi(theta_samples, hdi_prob=0.95)
array([0.39564857, 0.66586659])
Nel caso presente, l’intervallo di credibilità che lascia la stessa probabilità nelle due code e le regioni HPD sono quasi identici.
30.4. Principali distribuzioni coniugate#
Esistono altre combinazioni di verosimiglianza e distribuzione a priori le quali producono una distribuzione a posteriori che ha la stessa densità della distribuzione a priori. Sono elencate qui sotto le più note coniugazioni tra modelli statistici e distribuzioni a priori.
Per il modello Normale-Normale
, la distribizione iniziale è e la distribuzione finale è .Per il modello Poisson-gamma
, la distribizione iniziale è e la distribuzione finale è .Per il modello esponenziale
, la distribizione iniziale è e la distribuzione finale è .Per il modello uniforme-Pareto
, la distribizione iniziale è e la distribuzione finale è .
30.5. Commenti e considerazioni finali#
Nel Capitolo è stato spiegato come combinare le conoscenze a priori e le evidenze fornite dai dati per ottenere una stima della distribuzione di probabilità a posteriori del parametro
30.6. Watermark#
%load_ext watermark
%watermark -n -u -v -iv
Last updated: Sat Jun 17 2023
Python implementation: CPython
Python version : 3.11.3
IPython version : 8.12.0
matplotlib: 3.7.1
pymc : 5.5.0
arviz : 0.15.1
seaborn : 0.12.2
numpy : 1.24.3
pandas : 1.5.3
scipy : 1.10.1