Un’introduzione a Pymc e il linguaggio per descrivere modelli statistici

Introduzione a Pymc e il linguaggio per modelli statistici

Foto di Paul Nicklen su Treehugger

Nel nostro articolo precedente su perché la maggior parte degli esempi di inferenza bayesiana rappresenta in modo errato ciò che è, abbiamo chiarito un comune fraintendimento tra i principianti di statistica bayesiana. Vale a dire, il campo della statistica bayesiana NON È definito dal suo uso del teorema di Bayes, ma piuttosto dal suo uso di distribuzioni di probabilità per caratterizzare l’incertezza e considerare l’intera gamma di possibili risultati. Quindi, ad esempio, anziché venire informati che un determinato dispositivo medico è efficace al 95% nel rilevare una malattia dato il fatto che si sa con certezza assoluta di essere già infetti (che possiamo altrimenti chiamare il nostro Tasso di Veri Positivi (TPR)), possiamo considerare scenari in cui quel dispositivo è efficace solo al 24%, al 69% o al 91% nel rilevare la malattia e come quei numeri cambiano la probabilità che ci infettiamo dato che abbiamo dato positivo al test del dispositivo. Rifiutando di accettare stime puntuali (come il TPR statico) nei nostri modelli e favorendo invece un approccio probabilistico mediante la considerazione di un’ampia gamma di risultati, rimaniamo molto più fedeli alla visione bayesiana della probabilità.

Detto questo, la prossima fase del nostro percorso di apprendimento verso il raggiungimento di un pensiero di tipo bayesiano riguardo alla teoria della probabilità consiste nella costruzione di modelli statistici bayesiani che integrino il loro approccio probabilistico per caratterizzare l’incertezza. Uno degli strumenti migliori per realizzare questo obiettivo è il pacchetto Python PyMC, che è stato costruito appositamente per condurre inferenza bayesiana e costruire modelli di apprendimento automatico probabilistici.

Tuttavia, prima di iniziare a costruire modelli bayesiani, dobbiamo fare una breve deviazione e introdurre il linguaggio (cioè le notazioni matematiche) usato dagli esperti per formulare questi modelli. Capire il linguaggio della costruzione dei modelli è pratico per due motivi:

Motivo 1. Il primo motivo è che imparare come gli esperti, come scienziati e ricercatori, formulano modelli utilizzando la notazione dei modelli statistici applicati è direttamente correlato al modo in cui scriviamo codice eseguibile su PyMC. Utilizzare le notazioni dei modelli consente al modelizzatore statistico di comunicare le loro supposizioni sottostanti su un modello utilizzando poche linee di notazione. Sosteniamo che le notazioni siano una forma di comunicazione molto più semplice rispetto all’alternativa, che impone all’interlocutore di ricordare parole complicate come omoschedasticità per comprendere le caratteristiche del modello con cui stanno lavorando.

Motivo 2. L’altro motivo è che questo linguaggio si estende a tutti i domini, dall’arte alla storia, all’astronomia e alla biologia della conservazione (McElreath, 2020). Pertanto, essendo a proprio agio con il linguaggio della descrizione dei modelli, espanderemo inevitabilmente i confini entro i quali possiamo comprendere e comunicare conoscenze all’interno della comunità scientifica (McElreath, 2020).

Come introduzione generale alla notazione dei modelli statistici applicati, il prof. McElreath (2020) descrive alcuni principi come base per descrivere e comprendere modelli nel linguaggio colloquiale scientifico:

  1. Gli esperti riconoscono un insieme di variabili con cui lavorare, che talvolta sono osservabili. Queste vengono chiamate dati.a. Fenomeni non osservabili che derivano dai dati, come le medie, sono definiti come parametri.
  2. Gli esperti definiscono ogni variabile in base ad altre variabili O in base a una distribuzione di probabilità.
  3. La combinazione di variabili e delle loro distribuzioni di probabilità definisce un modello generativo congiunto, che possiamo utilizzare per simulare un’osservazione ipotetica e analizzare quelle reali.

Ora che abbiamo la ragione dietro il mondo delle notazioni dei modelli, illustreremo come le utilizziamo con un esempio giocattolo. Supponiamo che tu sia un conservazionista interessato a quantificare la probabilità di incontrare un lupo maschio della costa dell’Isola di Vancouver (Canis lupus crassodon) rispetto a una femmina in natura.

Come nota a margine, uno dei motivi per cui siamo interessati a utilizzare il lupo della costa dell’Isola di Vancouver nel nostro esempio è perché sono stati scoperti e classificati di recente come una sottospecie separata. Pertanto, la quantità di dati e informazioni su di loro è limitata rispetto ad altre sottospecie di lupi, quindi sarebbe interessante generare stime sulle caratteristiche di questa sottopopolazione qui e in futuri articoli. Non solo i dati sono limitati, ma a nostro parere, sono anche animali molto interessanti che si distinguono da altre sottospecie di lupi per la loro capacità di nuotare lunghe distanze e incorporare cibo di mare nella loro dieta naturale (Muñoz‐Fuentes et al., 2009). Da una prospettiva più ampia, la ricerca sull’ecologia animale e le statistiche di popolazione possono avere implicazioni sullo stato di conservazione di un animale, sulla gestione delle popolazioni e sulla comprensione del suo ruolo all’interno di un ecosistema più ampio.

Tornando al nostro esempio, basandoci solo su questa semplice domanda di ricerca, possiamo ora generare delle ipotesi e costruire un modello che quantifichi il rapporto di genere tra i lupi costieri dell’Isola di Vancouver. Per cominciare, possiamo assumere che la distribuzione dei risultati in questo esperimento mentale sia binomiale, il che significa che possono verificarsi solo due esiti come risultato della nostra domanda di ricerca. O ci imbattiamo in un lupo costiero maschio lungo il nostro viaggio attraverso l’Isola di Vancouver o in una femmina. Infatti, possiamo effettivamente simulare come dovrebbe apparire una distribuzione binomiale quando c’è una probabilità del 50% che si verifichi uno dei due esiti binomiali utilizzando il seguente codice:

# Importiamo le librerie necessarieimport arviz as azimport matplotlib.pyplot as pltimport numpy as npimport pandas as pdimport pymc as pmimport scipy.stats as statsaz.style.use("arviz-darkgrid")

x_binom = np.random.binomial(n=10, p=0.5, size=100)sns.histplot(x_binom, kde=True, discrete=True)plt.xlabel(“Numero di volte in cui la moneta è caduta su testa (N=10)”)plt.ylabel(“Numero di campioni”)plt.title('Visualizzazione di una distribuzione binomiale')

Dal grafico, otterremo un’approssimazione vicina a una distribuzione normale dal risultato di ciò che equivale a 10 lanci di moneta, con l’asse x che descrive il numero di volte in cui la moneta è caduta su testa. E se ripetiamo questo esperimento 100 volte come abbiamo fatto nel codice, l’asse y indica quante volte otteniamo il risultato sull’asse x. In altre parole, su 100 volte in cui abbiamo ripetuto questo esperimento, circa 30 di questi esperimenti hanno portato alla moneta che cade su testa circa 5 volte su 10 lanci. Nota che i risultati del grafico varieranno leggermente ogni volta che esegui il codice sopra.

Aggiungiamo anche un’altra ipotesi nel nostro modello, in cui assumiamo che entrambi gli esiti abbiano uguale probabilità di verificarsi. Tra le distribuzioni binomiali, questo potrebbe non essere sempre il caso, quindi non è un’ipotesi che dovremmo sempre prendere come predefinita. Le distribuzioni in cui c’è una possibilità uguale di qualsiasi risultato sono chiamate distribuzione uniforme. Viene spesso utilizzata nelle notazioni dei modelli come un modo per comunicare l’assenza di conoscenze precedenti sugli esiti possibili nel nostro modello. Tuttavia, con sufficiente pratica ed esperienza nel tempo, scopriremo che è spesso raro che un praticante bayesiano che costruisce il modello abbia conoscenze precedenti nulle sul fenomeno che intendono modellare. Possiamo anche simulare questa distribuzione se sei curioso di vedere come appare in un grafico:

x_uni = np.linspace(-.25, 1.25, 100)plt.plot(x_uni, stats.uniform.pdf(x_uni, 0, 1))plt.xlabel("Prior (Rapporto di genere)")plt.ylabel("Plausibilità del rapporto")plt.title("Visualizzazione di una distribuzione uniforme")

Una volta identificate le nostre ipotesi sul modello di rapporto di genere dei lupi costieri, formalizziamole nella notazione del modello statistico applicato:

  • M ~ Binomial(N, p): Possiamo interpretare la nostra prima riga come “il numero di maschi (M) che incontreremo è distribuito in modo binomiale con una dimensione del campione di N e una probabilità precedente di p.” Questa riga rappresenterà la nostra funzione di verosimiglianza, anche nota come probabilità della prova P(E | H) nel teorema di Bayes, che è la probabilità di osservare un dato risultato. In questo caso, il “dato risultato” sarà il risultato dei nostri dati campionati sui lupi costieri, di cui parleremo più dettagliatamente tra poco.
  • p ~ Uniform(0, 1): Infine, possiamo interpretare la nostra seconda riga come “la probabilità precedente P(H) è distribuita in modo uniforme tra 0 e 1.”

Si noti che il simbolo ~ significa “è distribuito come” e comunica che il parametro sul lato sinistro del modello ha una relazione stocastica (che significa “determinato casualmente” o ha una proprietà di “distribuzione di probabilità casuale” che può essere prevista statisticamente MA NON precisamente) con la sua distribuzione. In altre parole, quando vediamo il simbolo ~ che rappresenta una relazione stocastica, ci aspettiamo che l’output dell’affermazione matematica a destra sia una distribuzione catturata dalla variabile a sinistra (nel nostro caso, sia p che M). Esprimere relazioni stocastiche è in netto contrasto con le affermazioni matematiche che siamo abituati a vedere con un singolo numero (piuttosto che una distribuzione) come output di un’equazione che consiste in un simbolo =, come ad esempio una comune che vediamo spesso per le Regressioni Lineari, che è y = mx + b.

L’ultima cosa di cui abbiamo bisogno prima di poter tradurre il nostro modello in codice è generare un campione di osservazioni per la nostra funzione di verosimiglianza ( P(E | H) ), che rappresenta la probabilità di osservare il dato risultato degli eventi. Ad esempio, se utilizzassimo i dati di campione acquisiti da Muñoz-Fuentes e i suoi colleghi (2009) nello studio sui lupi costieri dell’Isola di Vancouver, scopriremmo che su 28 lupi di sesso identificato, 15 erano maschi. In questo caso, possiamo descrivere la nostra funzione di verosimiglianza come la probabilità di incontrare 15 maschi su un totale di 28 lupi che abbiamo incontrato nel nostro dataset di campione.

Con le nostre variabili stabilite, possiamo ora utilizzare PyMC per incorporare direttamente nel codice le notazioni del modello che abbiamo appena formalizzato:

# Leggendo il nostro dataset di campione dallo studio di Muñoz-Fuentes et al. (2009).wolf_samples = "https://raw.githubusercontent.com/vanislekahuna/Statistical-Rethinking-PyMC/main/Data/Vancouver-Island-wolf-samples.csv"islandwolves_df = pd.read_csv(wolf_samples)islandwolves_df# Determinazione del rapporto tra maschi e femmine nel datasetprint(f"Il conteggio dei valori per il sesso dei lupi costieri dell'Isola di Vancouver è: ")print(islandwolves_df["Sexb"].value_counts()

# Generazione delle osservazioni di campione per il nostro modello.males = islandwolves_df["Sexb"].value_counts()["M"]females = islandwolves_df["Sexb"].value_counts()["F"]total_samples = males + females

### Il nostro modello PyMC ###with pm.Model() as wolf_ratio:  # Una distribuzione uniforme di possibili risultati per i valori precedenti. priors = pm.Uniform("uniform_prior", lower=0, upper=1)   # Genera una distribuzione binomiale dalla nostra distribuzione uniforme dei valori precedenti e dai nostri dati di campione, in cui 15 dei 28 lupi si sono rivelati maschi.  n_males = pm.Binomial("sex_ratio", p=priors, observed=males, n=total_samples)  # Estrae campioni dal modello per generare la nostra distribuzione posteriore.   trace_wolves = pm.sample(1000, tune=1000)######################

Ora che il nostro modello è stato attualizzato, non perdiamoci nei dettagli e dimentichiamo lo scopo principale di costruire il modello in primo luogo. Il motivo per cui stiamo costruendo il modello è quello di quantificare la probabilità di incontrare un lupo maschio attraverso la distribuzione delle possibilità e di classificarli essenzialmente in base alla loro verosimiglianza. In altre parole, siamo interessati alla probabilità di osservare 15 maschi su un totale di 28 lupi che abbiamo incontrato (la nostra funzione di verosimiglianza) dato ogni possibile rapporto tra maschi e femmine, come rappresentato dalla nostra distribuzione precedente (ad esempio, il 1% di tutti i lupi selvatici sono maschi, il 2%, il 3%, fino al 100%). Il risultato finale qui dovrebbe essere una distribuzione di probabilità posteriore che mostra quali rapporti tra maschi e femmine sono più plausibili in base ai dati di campione che abbiamo osservato. Dal codice sopra, la distribuzione posteriore è catturata nell’oggetto `trace_wolves`, che genera le probabilità posteriori risultanti campionando il modello che abbiamo costruito. Nel caso ti stia chiedendo, i metodi di campionamento utilizzati per generare la distribuzione posteriore sono un argomento che possiamo trattare in un articolo successivo, poiché PyMC fornisce diverse opzioni, ognuna delle quali è un proprio profondo ed esteso labirinto.

Il grafico a sinistra di seguito visualizza come appare la nostra distribuzione posteriore risultante, con l’asse x che rappresenta la nostra distribuzione precedente di probabilità dei rapporti di genere e l’asse y che rappresenta la plausibilità di quel rapporto:

az.plot_trace(trace_wolves)

E se può essere utile strutturare il modello `wolf_ratio` che abbiamo appena costruito come variabili nel Teorema di Bayes, possiamo farlo:

Di nuovo, è importante sottolineare che l’obiettivo qui è generare una distribuzione posteriore di plausibilità per ogni valore precedente nella sua distribuzione. NON STIAMO cercando di trovare il rapporto di genere con la probabilità più alta. Detto questo, se fossi interessato a quale sarebbe quel valore, potremmo utilizzare la funzione `find_map` di PyMC per fornirci il valore precedente più plausibile (cioè il rapporto maschi/femmine più plausibile), nella distribuzione posteriore:

pm.find_MAP(model=wolf_ratio)

O la funzione `summary` di Arviz, che ci fornisce anche il massimo a priori (MAP) della distribuzione posteriore, oltre agli intervalli di densità più elevata (a volte chiamati anche intervalli di compatibilità), che sono i due valori precedenti che contengono l’89% della distribuzione.

az.summary(trace_wolves, round_to=2, kind="stats", hdi_prob=0.89)

Ecco fatto. Abbiamo appena fatto una breve introduzione al linguaggio delle notazioni dei modelli statistici utilizzati nella maggior parte degli studi accademici. Abbiamo quindi integrato le notazioni dei modelli nel nostro primo modello PyMC, in cui abbiamo simulato la plausibilità di incontrare un lupo costiero maschio dell’Isola di Vancouver in natura, dato che non sappiamo nulla sul rapporto tra maschi e femmine. Familiarizzarsi con le notazioni dei modelli formerà una base solida per la costruzione di modelli bayesiani più complessi utilizzando PyMC. Stiamo pianificando di pubblicare ulteriori articoli su questo argomento in futuro, a partire da un articolo sulla previsione lineare bayesiana, quindi rimanete sintonizzati!

Riferimenti:

McElreath, R. (2020). Statistical Rethinking: A Bayesian Course with examples in R and Stan. Routledge.

Muñoz‐Fuentes, V., Darimont, C. T., Wayne, R. K., Paquet, P. C., & Leonard, J. A. (2009). Ecological factors drive differentiation in wolves from British Columbia. Journal of Biogeography, 36(8), 1516–1531. https://www.raincoast.org/files/publications/papers/Munoz_et_al_2009_JBiogeog.pdf

Muñoz-Fuentes, V., Darimont, C.T., Paquet, P.C., & Leonard, J.A. (2010). The genetic legacy of extirpation and re-colonization in Vancouver Island wolves. Conservation Genetics, 11, 547–556. https://digital.csic.es/bitstream/10261/58619/1/evol.pdf