Modello di Poisson

Modello di Poisson#

In questo capitolo, esamineremo la stima del parametro “rate” di un modello di Poisson utilizzando PyMC. Il modello di Poisson è uno strumento impiegato per modellare il conteggio di eventi rari in un intervallo di tempo o spazio. Il parametro “rate” in questo contesto rappresenta il tasso di incidenza medio degli eventi nel sistema considerato, ovvero il numero medio di eventi che si verificano nell’unità di tempo o spazio.

Attraverso l’utilizzo di PyMC, otterremo la distribuzione a posteriori del parametro “rate” che tiene conto dell’incertezza nella stima. Ciò ci permetterà di ottenere stime puntuali e intervalli di credibilità per il parametro “rate”, fornendo una migliore comprensione dell’andamento degli eventi nel sistema considerato e una misura dell’incertezza associata alla nostra stima.

Il vantaggio principale di utilizzare la stima MCMC con PyMC è la sua capacità di gestire casi in cui i dati possono essere scarsi o non soddisfare gli assunti delle tecniche di stima tradizionali. Inoltre, attraverso questo approccio, potremo esplorare il comportamento del modello di Poisson e valutarne la capacità di descrivere correttamente il fenomeno oggetto di studio, consentendoci di ottenere una migliore comprensione e interpretazione dei dati osservati.

Preparazione del Notebook#

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import scipy.stats as stats
import pymc as pm
import pymc.sampling_jax
import arviz as az
import warnings

warnings.filterwarnings("ignore", category=UserWarning)
warnings.filterwarnings("ignore", category=FutureWarning)
warnings.filterwarnings("ignore", category=Warning)

/Users/corrado/opt/anaconda3/envs/pymc_env/lib/python3.11/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
  from .autonotebook import tqdm as notebook_tqdm

%config InlineBackend.figure_format = 'retina'
RANDOM_SEED = 42
rng = np.random.default_rng(RANDOM_SEED)
az.style.use("arviz-darkgrid")
sns.set_theme(palette="colorblind")

Domanda della ricerca#

Come spiegato qui, i dati che esamineremo sono raccolti dal Washington Post con lo scopo di registrare ogni sparatoria mortale negli Stati Uniti ad opera di agenti di polizia, a partire dal 1° gennaio 2015. Il Washington Post ha adottato un approccio sistematico e accurato nella raccolta di queste informazioni, fornendo dati che possono essere utili per valutare i problemi legati alla violenza delle forze di polizia negli Stati Uniti.

Lo scopo della presente analisi dei dati è determinare il tasso di sparatorie fatali da parte della polizia negli Stati Uniti per ogni anno, e fornire una stima dell’incertezza associata a questo valore.

Importazione e pre-processing dei dati#

url = "https://raw.githubusercontent.com/washingtonpost/data-police-shootings/master/v2/fatal-police-shootings-data.csv"
fps_dat = pd.read_csv(url)
fps_dat.head()

Show code cell output Hide code cell output

	id	date	threat_type	flee_status	armed_with	city	county	state	latitude	longitude	location_precision	name	age	gender	race	race_source	was_mental_illness_related	body_camera	agency_ids
0	3	2015-01-02	point	not	gun	Shelton	Mason	WA	47.246826	-123.121592	not_available	Tim Elliot	53.0	male	A	not_available	True	False	73
1	4	2015-01-02	point	not	gun	Aloha	Washington	OR	45.487421	-122.891696	not_available	Lewis Lee Lembke	47.0	male	W	not_available	False	False	70
2	5	2015-01-03	move	not	unarmed	Wichita	Sedgwick	KS	37.694766	-97.280554	not_available	John Paul Quintero	23.0	male	H	not_available	False	False	238
3	8	2015-01-04	point	not	replica	San Francisco	San Francisco	CA	37.762910	-122.422001	not_available	Matthew Hoffman	32.0	male	W	not_available	True	False	196
4	9	2015-01-04	point	not	other	Evans	Weld	CO	40.383937	-104.692261	not_available	Michael Rodriguez	39.0	male	H	not_available	False	False	473

# Convert date
fps_dat["date"] = pd.to_datetime(fps_dat["date"])

# Create a new column 'year' to store the year information from the 'date' column
fps_dat["year"] = fps_dat["date"].dt.year

fps_dat.columns

Index(['id', 'date', 'threat_type', 'flee_status', 'armed_with', 'city',
       'county', 'state', 'latitude', 'longitude', 'location_precision',
       'name', 'age', 'gender', 'race', 'race_source',
       'was_mental_illness_related', 'body_camera', 'agency_ids', 'year'],
      dtype='object')

# Filter out rows with year equal to 2023
fps = fps_dat[fps_dat["year"] != 2023]

# Count occurrences of each year in fps
year_counts = fps["year"].value_counts()
print(year_counts)

year
2022    1096
2021    1048
2020    1019
2019     997
2015     995
2018     992
2017     983
2016     958
2024      60
Name: count, dtype: int64

# Define the data
data = {
    "year": [2022, 2021, 2020, 2019, 2015, 2018, 2017, 2016],
    "events": [1096, 1048, 1019, 997, 995, 992, 983, 958],
}
# Convert the data to a DataFrame
df = pd.DataFrame(data)
print(df)

   year  events
0  2022    1096
1  2021    1048
2  2020    1019
3  2019     997
4  2015     995
5  2018     992
6  2017     983
7  2016     958

Modello di Poisson#

Il nostro interesse riguarda il tasso di occorrenza di sparatorie fatali da parte della polizia per anno. Indicheremo questo tasso come $θ$ , e il suo intervallo di valori possibili è $[0, \infty)$ . Un modello di Poisson rappresenta tipicamente il punto di partenza per l’analisi di dati relativi alle frequenze assolute di un evento in un intervallo di tempo fissato. Il modello presuppone che i dati seguano una distribuzione di Poisson con un parametro di tasso $θ$ :

P (y ∣ θ) \propto θ^{y} * \exp (- θ),

dove i dati $y$ possono assumere solo valori interi non negativi.

Distribuzione a priori#

Come distribuzione a priori per il parametro $θ$ nel modello di Poisson possiamo usare la distribuzione Gamma, poiché è una scelta coniugata. Ciò significa che, quando viene combinata con la distribuzione di Poisson come verosimiglianza dei dati, la distribuzione Gamma produce una distribuzione a posteriori con una forma analitica semplice. Questa caratteristica semplifica il processo di inferenza bayesiana.

Nel nostro caso, il parametro $θ$ rappresenta il tasso di occorrenza di sparatorie fatali per anno negli Stati Uniti. Prima di osservare i dati effettivi riportati dal Washington Post, abbiamo una conoscenza limitata su tale fenomeno. Pertanto, dobbiamo specificare una distribuzione a priori per $θ$ che rifletta la nostra incertezza iniziale. As esempio, possiamo ipotizzare che ci sia, in media, una sparatoria mortale per stato al mese, quindi 12 sparatorie all’anno per stato. Questo ci porta a una stima iniziale di 600 sparatorie fatali negli Stati Uniti ogni anno. Dato che non siamo molto sicuri di questa ipotesi, vogliamo specificare una distribuzione a priori con un certo grado di incertezza. Imponiamo dunque una deviazione standard pari a 200.

Creiamo un istogramma della distribuzione Gamma con i parametri specificati sopra usando PyMC.

mu = 600  # mu hyperparameter for the Gamma prior
sigma = 200  # sigma hyperparameter for the Gamma prior

x = pm.Gamma.dist(mu=mu, sigma=sigma)
x_draws = pm.draw(x, draws=50000, random_seed=2)
sns.histplot(x_draws)

plt.xlabel("Tasso di occorrenza (anno)")
plt.ylabel("Frequenza")
plt.title("Distribuzione Gamma con mu = 600 e sigma = 200")
plt.show()

Modello di Poisson con PyMC#

In PyMC, la distribuzione Gamma è parametrizzata con alpha = shape e beta = 1/scale. Nel nostro caso, shape = 3 e scale = 200. Pertanto, nella parametrizzazione PyMC, la distribuzione Gamma avrà parametri alpha = 3 e beta = 1 / 200. Formuliamo dunque il modello di Poisson usando questi iper-parametri per la distribuzione a priori del parametro $λ$ (rate) della distribuzione di Poisson.

with pm.Model() as poisson_model:
    rate = pm.Gamma("rate", mu=mu, sigma=sigma)
    # Likelihood (Poisson) using the observed data
    y = df["events"].values
    events = pm.Poisson("events", mu=rate, observed=y)

Eseguiamo il prior predictive check per verificare che l’instanziazione del modello sia corretta per quel che riguarda la distribuzione a priori.

# Use your model to sample from the prior
with poisson_model:
    idata_prior = pm.sample_prior_predictive()

az.plot_trace(idata_prior.prior);

La distribuzione predittiva a priori suggerisce che la distribuzione a priori per il parametro rate è adeguata per i dati a disposizione. Eseguiamo dunque il campionamento.

with poisson_model:
    idata = pm.sampling_jax.sample_numpyro_nuts()

Show code cell output Hide code cell output

Compiling...

Compilation time = 0:00:01.096445

Sampling...

  0%|                                                                                  | 0/2000 [00:00<?, ?it/s]

Compiling.. :   0%|                                                                    | 0/2000 [00:00<?, ?it/s]

  0%|                                                                                  | 0/2000 [00:00<?, ?it/s]

Compiling.. :   0%|                                                                    | 0/2000 [00:00<?, ?it/s]

  0%|                                                                                  | 0/2000 [00:00<?, ?it/s]

Compiling.. :   0%|                                                                    | 0/2000 [00:00<?, ?it/s]

  0%|                                                                                  | 0/2000 [00:00<?, ?it/s]

Compiling.. :   0%|                                                                    | 0/2000 [00:00<?, ?it/s]

Running chain 1:   0%|                                                                 | 0/2000 [00:01<?, ?it/s]

Running chain 3:   0%|                                                                 | 0/2000 [00:01<?, ?it/s]

Running chain 0:   0%|                                                                 | 0/2000 [00:01<?, ?it/s]

Running chain 2:   0%|                                                                 | 0/2000 [00:01<?, ?it/s]

Running chain 0: 100%|████████████████████████████████████████████████████| 2000/2000 [00:01<00:00, 1932.50it/s]

Running chain 1: 100%|████████████████████████████████████████████████████| 2000/2000 [00:01<00:00, 1934.62it/s]

Running chain 2: 100%|████████████████████████████████████████████████████| 2000/2000 [00:01<00:00, 1937.39it/s]

Running chain 3: 100%|████████████████████████████████████████████████████| 2000/2000 [00:01<00:00, 1940.15it/s]

Sampling time = 0:00:01.274469

Transforming variables...

Transformation time = 0:00:00.066985

Esaminiamo la distribuzione a posteriori di rate.

az.plot_trace(idata, combined=True, kind="rank_bars");

Il modello converge rapidamente e i grafici delle tracce sembrano ben mescolati.

Generiamo un sommario numerico della distribuzione a posteriori.

pm.summary(idata)

	mean	sd	hdi_3%	hdi_97%	mcse_mean	mcse_sd	ess_bulk	ess_tail	r_hat
rate	1010.358	10.908	989.279	1030.241	0.276	0.196	1555.0	2085.0	1.0

Usiamo ArviZ per generare l’intervallo di credibilità al 94% per la distribuzione a posteriori del parametro rate.

az.plot_posterior(idata, var_names="rate");

Eseguiamo il posterior predictive check.

with poisson_model:
    pm.sample_posterior_predictive(idata, extend_inferencedata=True)

Sampling: [events]

100.00% [4000/4000 00:00<00:00]

Possiamo utilizzare az.plot_ppc() per verificare che i campioni predittivi a posteriori siano simili ai dati osservati.

az.plot_ppc(idata, num_pp_samples=300);

In conclusione, possiamo affermare che il tasso stimato di sparatorie fatali da parte della polizia negli Stati Uniti è di 1011 casi all’anno, con un intervallo di credibilità del 94% di [988.53, 1031.41].

Derivazione analitica#

In realtà, il problema precedente poteva essere risolto in maniera analitica, senza ricorrere al campionamento MCMC. Infatti, la distribuzione coniugata al modello di Poisson è una densità di tipo gamma e dunque la distribuzione a posteriori è ancora una distribuzione gamma.

Consideriamo un modello di Poisson con verosimiglianza:

P (Y | θ) = \frac{e^{- θ} \cdot θ^{Y}}{Y!},

dove $Y$ è il numero di eventi osservati e $θ$ è il parametro che rappresenta il tasso di incidenza degli eventi. Supponiamo di avere una distribuzione a priori gamma con i parametri $a$ e $b$ :

P (θ) = \frac{b^{a}}{Γ (a)} \cdot θ^{a - 1} \cdot e^{- b θ},

dove $a$ e $b$ sono i parametri della distribuzione gamma. Per trovare la distribuzione a posteriori, dobbiamo moltiplicare la verosimiglianza e la distribuzione a priori, e normalizzare il risultato in modo che la distribuzione risultante sia una distribuzione di probabilità valida. La distribuzione a posteriori è quindi data da:

P (θ | Y) \propto P (Y | θ) \cdot P (θ) = \frac{e^{- θ} \cdot θ^{Y}}{Y!} \cdot \frac{b^{a}}{Γ (a)} \cdot θ^{a - 1} \cdot e^{- b θ}

Per ottenere la distribuzione a posteriori normalizzata, dobbiamo dividere il risultato ottenuto per l’integrale di questa espressione su tutto l’intervallo dei valori di $θ$ . Semplificando e normalizzando questa espressione, otteniamo la distribuzione a posteriori per $θ$ uguale ad una distribuzione gamma con parametri aggiornati, ovvero:

P (θ | Y) = \frac{b^{' a^{'}}}{Γ (a^{'})} \cdot θ^{a^{'} - 1} \cdot e^{- b^{'} θ},

dove i parametri aggiornati sono dati da

\begin{array}{r} a^{'} = a + Y \\ b^{'} = b + 1 \end{array}

In conclusione, la distribuzione a posteriori per $θ$ , data una verosimiglianza di Poisson e una distribuzione a priori gamma, è ancora una distribuzione gamma con parametri aggiornati.

Risolviamo dunque il problema precedente utilizzando la formula della distribuzione a posteriori riportata sopra. Se la distribuzione a priori è una distribuzione gamma di parametri α e β, la distribuzione a posteriori è ancora una gamma con parametri aggiornati. Possiamo calcolare i parametri della distribuzione a posteriori nel modo seguente:

Parametro di forma a posteriori (α_post) = α_prior + Σ(y_i), dove Σ(y_i) rappresenta la somma dei dati osservati.
Parametro di tasso a posteriori (β_post) = β_prior + n, dove n è il numero di punti dati.

Con questi parametri aggiornati, possiamo poi calcolare la media a posteriori della distribuzione gamma e l’intervallo di credibilità.

data = [1096, 1048, 1019, 997, 995, 992, 983, 958]

# Prior hyperparameters
alpha_prior = 3
beta_prior = 1 / 200

# Data summary
n = len(data)
sum_y = np.sum(data)

# Posterior hyperparameters
alpha_post = alpha_prior + sum_y
beta_post = beta_prior + n

# Posterior distribution (Gamma)
posterior_gamma = stats.gamma(alpha_post, scale=1 / beta_post)

# Calculate the mean and credibility interval (94%)
posterior_mean = posterior_gamma.mean()
credible_interval = posterior_gamma.interval(0.94)

print("Estimated Rate (Posterior Mean):", posterior_mean)
print("Credibility Interval (94%):", credible_interval)

Estimated Rate (Posterior Mean): 1010.7432854465958
Credibility Interval (94%): (989.7152329772036, 1031.9826527942685)

L’output delle istruzioni precedenti fornisce il tasso stimato a posteriori e l’intervallo di credibilità al 94%. A causa di approssimazioni numeriche, i valori non coincidono esattamente con i risultati ottenuti con PyMC, ma sono molto simili.

Un altro esempio di applicazione del modello di Poisson è presentato in appendice.

%run ../wm.py

Watermark:
----------
Last updated: 2024-01-26T18:58:51.293717+01:00

Python implementation: CPython
Python version       : 3.11.7
IPython version      : 8.19.0

Compiler    : Clang 16.0.6 
OS          : Darwin
Release     : 23.3.0
Machine     : x86_64
Processor   : i386
CPU cores   : 8
Architecture: 64bit

Modello di Poisson

Contents

Modello di Poisson#

Preparazione del Notebook#

Domanda della ricerca#

Importazione e pre-processing dei dati#

Modello di Poisson#

Distribuzione a priori#

Modello di Poisson con PyMC#

Derivazione analitica#