La predizione delle frequenze

38. La predizione delle frequenze#

La predizione delle frequenze si riferisce al problema di prevedere il numero di volte che un evento si verificherà in futuro, sulla base delle informazioni disponibili. Questo problema è spesso affrontato utilizzando il teorema di Bayes, che consente di calcolare la probabilità di un evento sulla base delle informazioni precedenti e delle nuove informazioni acquisite. Per fare ciò, si utilizza una distribuzione a priori per rappresentare le nostre credenze iniziali sulla frequenza dell’evento e si aggiorna questa distribuzione utilizzando le nuove informazioni acquisite. La distribuzione risultante, ovvero la distribuzione a posteriori, rappresenta le nostre credenze aggiornate sulla frequenza dell’evento e può essere utilizzata per fare previsioni sulle future occorrenze dell’evento.

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import scipy as sc
import statistics
import arviz as az
import seaborn as sns
import scipy.stats as st
import pymc as pm
import xarray as xr
import warnings

warnings.filterwarnings("ignore")
warnings.simplefilter("ignore")

print(f"Running on PyMC v{pm.__version__}")

Running on PyMC v5.5.0

%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"

38.1. La ricerca sul trauma#

Per fare un esempio di questo approccio, in questo capitolo faremo riferimento alla ricerca sul trauma. In questo campo di ricerca, come in altri, il risultato di interesse può essere rappresentato dalla frequenza del numero di episodi che si verificano in un dato periodo di tempo. Nel caso della ricerca sulla violenza domestica, ad esempio, potremmo esaminare il tasso di atti aggressivi durante l’intervallo tra un momento temporale di baseline e un’intervista di follow-up. Altri esempi di risultati esprimibili in termini di frequenze nalla ricerca post-traumatica includono la frequenza dell’abuso di sostanze durante un periodo di osservazione o il numero di interventi della polizia durante un dato periodo.

Nell’esempio seguente esamineremo l’uso della predizione bayesiana per predire il numero di aggressioni nei confronti del partner nelle relazioni di coppia. I dati presentati sono tratti da uno studio che esamina la frequenza di episodi di aggressione messi in atto da pazienti di sesso maschile che avevano recentemente iniziato un programma di trattamento dell’alcol nei confronti del loro partner di sesso femminile.

La frequenza degli episodi di violenza, così come altri fenomeni quantificabili in termini di frequenze assolute, può essere modellata da un processo di Poisson, che si basa sul presupposto che gli eventi siano casuali e abbiano la stessa probabilità di verificarsi in qualsiasi momento. Naturalmente, questa ipotesi non è sempre valida, ma spesso è sufficiente per la modellazione.

38.2. La distribuzione a priori#

In una ricerca sui pazienti che avevano recentemente iniziato un programma di trattamento per l’abuso o la dipendenza da alcol, Gagnon et al. [GDLB+08] hanno trovato che, in un periodo di 6 mesi, il numero di assalti fisici da parte dei pazienti di genere maschile nei confronti del loro partner femminile è uguale, in media, a 11.46 (\(SD\) = 25.79; \(n\) = 114).

Per questi dati, è dunque sensato pensare che la distribuzione del numero di episodi di aggressione fisica può essere rappresentata dalla distribuzione esponenziale.

Dal punto di vista statistico, ricordiamo che la distribuzione esponenziale modella il numero di eventi che si verificano in un intervallo di tempo quando questi eventi si verificano raramente e indipendentemente l’uno dall’altro.

La funzione di densità di probabilità della distribuzione esponenziale è data da:

\[ f(x) = λe^{(-λx)} \]

dove λ è il parametro della distribuzione e rappresenta il tasso medio di occorrenza degli eventi. La media e la varianza della distribuzione esponenziale sono entrambe uguali a 1/λ.

Poniamoci dunque il problema di rappresentare le nostre credenze a priori relative al numero di episodi di aggressione fisica mediante una distribuzione esponenziale.

Dalla ricerca di Gagnon et al. [GDLB+08] sappiamo che, in un periodo di 6 mesi, il numero di assalti fisici nei confronti del partner femminile, in questa popolazione, è uguale a 11.46.

Considereremo una distribuzione esponenziale per rappresentare le nostre credenze a priori circa la frequenza media \(\mu\) degli episodi di aggressione nei confronti del partner per questa popolazione, in un periodo di 6 mesi. In Python, scipy.stats.expon è un modulo che fornisce funzioni per lavorare con la distribuzione esponenziale. In particolare, la funzione pdf (probability density function) calcola la funzione di densità di probabilità della distribuzione esponenziale per un dato valore di x.

La sintassi per utilizzare questa funzione è la seguente:

scipy.stats.expon.pdf(x, loc=0, scale=1)

dove x è il valore per il quale si vuole calcolare la funzione di densità di probabilità. Il parametro loc(opzionale) e scale specificano rispettivamente la posizione e la scala della distribuzione. La posizione (loc) di solito è impostata a 0, mentre la scala (scale) è l’inverso del parametro λ della distribuzione esponenziale. In altre parole, la scala è uguale alla media della distribuzione esponenziale.

Per il caso presente, dunque, se vogliamo che la distribuzione esponenziale abbia una media di 11.46, procediamo come indicato sotto

x = np.linspace (0, 50, 100) 
y = st.expon.pdf(x, 0, 11.46)
plt.plot(x, y)

[<matplotlib.lines.Line2D at 0x128e926d0>]

_images/0f1c1cb4f4cf0c5e73ff8de57bdd00425e3b0952504df358010718c4efbec490.svg

Per verificare che abbiamo implementato correttamente la funzione esponenziale con il parametro voluto, estraiamo un grande numero di realizzazioni della v.c. e calcoliamo la media.

r = st.expon.rvs(0, 11.46, size=100000)
r.mean()

11.45794424722882

38.3. Inferenza bayesiana#

Una volta capito come descrivere le nostre credenze a priori, poniamoci il problema di usare PyMC per l’inferenza Bayesiana.

Consideriamo un singlo individuo di genere maschile appartenente a questa popolazione. Se, in media, in 6 mesi ci aspettiamo un numero di episodi di violenza pari a 11.46, possiamo descrivere il numero di episodi di violenza per un singolo individuo con la seguente distribuzione esponenziale. Si noti che, in questo caso, la funzione pm.Exponential è parametrizzata usando il parametro l che è uguale a \(\lambda = 1/ \mu\).

l = 1/11.46

with pm.Model() as model:
    mu = pm.Exponential("mu", l)
    idata = pm.sample_prior_predictive(samples=10000, random_seed=rng)

Sampling: [mu]

Esaminiamo 10000 campioni casuali estratti dalla distribuzione a priori. Il risultato è simile alla distribuzione di densità teorica rappresentata nel grafico precedente.

_ = az.plot_posterior(idata.prior.mu)

_images/cfbf38c5cc72f0e0fd8a73923e606f4161171fb0481dc4b7d692f22722dacad5.svg

Si noti che, in questo caso, stiamo usando una distribuzione a priori informativa.

38.4. La verosimiglianza#

Adesso inseriamo nel modello PyMC la verosimiglianza.

In questo caso, usiamo quale modello generativo dei dati una distribuzione di Poisson.

Le distribuzioni di Poisson sono usate per modellare il numero di eventi rari che si verificano in un intervallo di tempo fisso. Ad esempio, il numero di episodi di comportamento inappropriato per settimana in individui con disturbi alimentari, nascite in un giorno o incidenti in una settimana.

La verosimiglianza di Poisson è simile a quella binomiale, ma non ha un limite superiore al numero di successi.

Per esempio, consideriamo un paziente chiamato Mario. Usando gli item della sottoscala relativa agli episodi di violenza fisica della Conflict Tactics Scales-2, troviamo che Mario ha avuto 8 episodi violenti nei confronti del partner negli ultimi 6 mesi.

Inseriamo questa informazione nel modello bayesiano usando 8 come dato che specifica una verosimiglianza di Poisson con parametro sconosciuto mu, a cui abbiamo attribuito una distribuzione a priori esponenziale ed eseguiamo il campionamento.

with pm.Model() as model:
    mu = pm.Exponential("mu", l)
    episodes = pm.Poisson("episodes", mu, observed=8)
    idata2 = pm.sample(2000)

Auto-assigning NUTS sampler...

Initializing NUTS using jitter+adapt_diag...

Multiprocess sampling (4 chains in 4 jobs)

NUTS: [mu]

Bad value in file PosixPath('/Users/corrado/mambaforge/envs/pymc_env/lib/python3.11/site-packages/arviz/plots/styles/arviz-whitegrid.mplstyle'), line 40 ("axes.prop_cycle: cycler('color', ['2a2eec', 'fa7c17', '328c06', 'c10c90', '933708', '65e5f3', 'e6e135', '1ccd6a', 'bd8ad5', 'b16b57'])"): Key axes.prop_cycle: "cycler('color', ['2a2eec', 'fa7c17', '328c06', 'c10c90', '933708', '65e5f3', 'e6e135', '1ccd6a', 'bd8ad5', 'b16b57'])" is not a valid cycler construction: 

---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
Cell In[7], line 4
mu = pm.Exponential("mu", l)
episodes = pm.Poisson("episodes", mu, observed=8)
----> 4 idata2 = pm.sample(2000)

File ~/mambaforge/envs/pymc_env/lib/python3.11/site-packages/pymc/sampling/mcmc.py:766, in sample(draws, tune, chains, cores, random_seed, progressbar, step, nuts_sampler, initvals, init, jitter_max_retries, n_init, trace, discard_tuned_samples, compute_convergence_checks, keep_warning_stat, return_inferencedata, idata_kwargs, nuts_sampler_kwargs, callback, mp_ctx, model, **kwargs)
_print_step_hierarchy(step)
try:
--> 766     _mp_sample(**sample_args, **parallel_args)
except pickle.PickleError:
   _log.warning("Could not pickle model, sampling singlethreaded.")

File ~/mambaforge/envs/pymc_env/lib/python3.11/site-packages/pymc/sampling/mcmc.py:1141, in _mp_sample(draws, tune, step, chains, cores, random_seed, start, progressbar, traces, model, callback, mp_ctx, **kwargs)
# We did draws += tune in pm.sample
draws -= tune
-> 1141 sampler = ps.ParallelSampler(
   draws=draws,
   tune=tune,
   chains=chains,
   cores=cores,
   seeds=random_seed,
   start_points=start,
   step_method=step,
   progressbar=progressbar,
   mp_ctx=mp_ctx,
)
try:
   try:

File ~/mambaforge/envs/pymc_env/lib/python3.11/site-packages/pymc/sampling/parallel.py:402, in ParallelSampler.__init__(self, draws, tune, chains, cores, seeds, start_points, step_method, progressbar, mp_ctx)
if mp_ctx.get_start_method() != "fork":
   step_method_pickled = cloudpickle.dumps(step_method, protocol=-1)
--> 402 self._samplers = [
   ProcessAdapter(
       draws,
       tune,
       step_method,
       step_method_pickled,
       chain,
       seed,
       start,
       mp_ctx,
   )
   for chain, seed, start in zip(range(chains), seeds, start_points)
]
self._inactive = self._samplers.copy()
self._finished: List[ProcessAdapter] = []

File ~/mambaforge/envs/pymc_env/lib/python3.11/site-packages/pymc/sampling/parallel.py:403, in <listcomp>(.0)
if mp_ctx.get_start_method() != "fork":
   step_method_pickled = cloudpickle.dumps(step_method, protocol=-1)
self._samplers = [
--> 403     ProcessAdapter(
       draws,
       tune,
       step_method,
       step_method_pickled,
       chain,
       seed,
       start,
       mp_ctx,
   )
   for chain, seed, start in zip(range(chains), seeds, start_points)
]
self._inactive = self._samplers.copy()
self._finished: List[ProcessAdapter] = []

File ~/mambaforge/envs/pymc_env/lib/python3.11/site-packages/pymc/sampling/parallel.py:259, in ProcessAdapter.__init__(self, draws, tune, step_method, step_method_pickled, chain, seed, start, mp_ctx)
   step_method_send = step_method
self._process = mp_ctx.Process(
   daemon=True,
   name=process_name,
   (...)
   ),
)
--> 259 self._process.start()
# Close the remote pipe, so that we get notified if the other
# end is closed.
remote_conn.close()

File ~/mambaforge/envs/pymc_env/lib/python3.11/multiprocessing/process.py:121, in BaseProcess.start(self)
assert not _current_process._config.get('daemon'), \
      'daemonic processes are not allowed to have children'
_cleanup()
--> 121 self._popen = self._Popen(self)
self._sentinel = self._popen.sentinel
# Avoid a refcycle if the target function holds an indirect
# reference to the process object (see bpo-30775)

File ~/mambaforge/envs/pymc_env/lib/python3.11/multiprocessing/context.py:300, in ForkServerProcess._Popen(process_obj)
@staticmethod
def _Popen(process_obj):
   from .popen_forkserver import Popen
--> 300     return Popen(process_obj)

File ~/mambaforge/envs/pymc_env/lib/python3.11/multiprocessing/popen_forkserver.py:35, in Popen.__init__(self, process_obj)
def __init__(self, process_obj):
   self._fds = []
---> 35     super().__init__(process_obj)

File ~/mambaforge/envs/pymc_env/lib/python3.11/multiprocessing/popen_fork.py:19, in Popen.__init__(self, process_obj)
self.returncode = None
self.finalizer = None
---> 19 self._launch(process_obj)

File ~/mambaforge/envs/pymc_env/lib/python3.11/multiprocessing/popen_forkserver.py:58, in Popen._launch(self, process_obj)
self.finalizer = util.Finalize(self, util.close_fds,
                              (_parent_w, self.sentinel))
with open(w, 'wb', closefd=True) as f:
---> 58     f.write(buf.getbuffer())
self.pid = forkserver.read_signed(self.sentinel)

KeyboardInterrupt: 

38.5. La distribuzione a posteriori#

Estraiamo dall’oggetto idata2 i campioni della distribuzione a posteriori del parametro mu (media del numero di episodi di violenza negli ultimi 6 mesi).

sample_posterior = idata2.posterior['mu']

Generiamo un grafico della distribuzione a posteriori del parametro mu.

az.plot_posterior(sample_posterior)

<Axes: title={'center': 'mu'}>

_images/926977290877f149fd48d1c4475b0b14a96d4116a6cc9399d09f885a890e8ce4.svg

Possiamo dunque concludere, con un livello di certezza soggettiva del 94%, che per Mario, il numero di episodi di violenza nei confronti del partner varierà da un minimo di 3.5 ad un massimo di 13, in un periodo di 6 mesi, con una media di 8.2.

Supponiamo ora di volere confrontare due individui, Mario e Paolo. Di Mario abbiamo osservato 8 episodi di violenza in 6 mesi; di Paolo abbiamo osservato 12 episodi di violenza negli ultimi 6 mesi. Possiamo dire che Paolo è più violento di Mario? Oppure dobbiamo pensare che la differenza tra i due sia solo una fluttuazione casuale?

Scriviamo il modello bayesiano nel modo seguente, usando sempre la distribuzione a priori che abbiamo definito in precedenza.

with pm.Model() as model3:
    mu_A = pm.Exponential("mu_A", l)
    mu_B = pm.Exponential("mu_B", l)
    episodes_A = pm.Poisson("episodes_A", mu_A, observed=[8])
    episodes_B = pm.Poisson("episodes_B", mu_B, observed=[12])
    idata3 = pm.sample(2000)

Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [mu_A, mu_B]

100.00% [12000/12000 00:02<00:00 Sampling 4 chains, 0 divergences]

Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 21 seconds.

Esaminiamo le distribuzioni a posteriori dei parametri mu_A e mu_B che rappresentano la media del numero di episodi di violenza per i due individui.

with model3:
    az.plot_trace(idata3, kind="rank_bars")
plt.tight_layout()

_images/ca95d91bcb0c2b205ed3c5aee06b0dec67036d7ccefcb9c4e2b5611a4eca70b0.svg

Estraiamo le distribuzioni a posteriori dei due parametri da idata3.

mu_A = idata3.posterior['mu_A']
mu_B = idata3.posterior['mu_B']
mu_B.mean(), mu_A.mean()

(<xarray.DataArray 'mu_B' ()>
 array(11.93233794),
 <xarray.DataArray 'mu_A' ()>
 array(8.21862358))

Rappresentiamo graficamente le due distribuzioni a posteriori.

az.plot_posterior(mu_A)
az.plot_posterior(mu_B)

<Axes: title={'center': 'mu_B'}>

_images/c9b236979314ed15d65828b2697a94628c156a245883796cce5977c5be417099.svg

_images/35b23246e70beeaa7f21be4535e565ae991bf0f98adffcba2256216f2bdcb477.svg

Vogliamo eseguire un test di ipotesi bayesiano per determinare la probabilità che la media del numero di episodi di violenza di Mario sia maggiore di quella di Paolo. Per fare ciò, calcoliamo quante volte mu_B è maggiore di mu_A nelle due distribuzioni a posteriori.

(mu_B > mu_A).mean()

<xarray.DataArray ()>
array(0.810875)

Possiamo dunque dire che, se confrontiamo i valori dei parametri delle due distribuzioni a posteriori, nell’81% di casi risulta che Paolo è più violento di Mario.

Il grafico seguente mostra le stime degli intervalli di credibilità del 94% per ciascuna delle 4 catene, pe i due parametri.

_ = az.plot_forest(idata3, var_names=["mu_A", "mu_B"])

_images/fefc9754a908a73493113f3ded6f85d9d3596b5cf507b41ed379d5506e9be56f.svg

Dato che gli intervalli di credibilità sono sovrapposti, concludiamo che non ci sono evidenze credibili di una differenza. Ovvero, sulla base delle nostre credenze a priori e sulla base dei dati osservati, ad un livello di certezza soggettiva del 94%, non possiamo concludere che Paolo sia più violento di Mario.

38.7. Watermark#

%load_ext watermark
%watermark -n -u -v -iv -w

Last updated: Sat May 06 2023

Python implementation: CPython
Python version       : 3.11.3
IPython version      : 8.13.2

numpy     : 1.23.5
xarray    : 2023.4.2
pandas    : 1.5.3
pymc      : 5.3.0
arviz     : 0.15.1
seaborn   : 0.12.2
scipy     : 1.10.1
matplotlib: 3.7.1

Watermark: 2.3.1