40. Confronto tra gruppi

L’obiettivo di questo capitolo è di ampliare la discussione del capitolo Inferenza bayesiana su una media, affrontando il confronto tra le medie di due o più gruppi indipendenti. Per cominciare, carichiamo le librerie necessarie.

from pymc import HalfNormal, HalfCauchy, Model, Normal, sample
import arviz as az
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pymc as pm
import seaborn as sns
import xarray as xr
import warnings

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

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

40.1. Stima della media di un gruppo

Iniziamo con l’utilizzo del caso precedentemente trattato nel capitolo precedente, in cui abbiamo calcolato l’incertezza relativa alla media di un singolo gruppo. Qui useremo i dati relativi ai pinguini Palmer, che verranno letti da un file csv, escludendo quelli che presentano valori mancanti:

penguins = pd.read_csv("data/penguins.csv")
# Subset to the columns needed
missing_data = penguins.isnull()[
    ["bill_length_mm", "flipper_length_mm", "sex", "body_mass_g"]
].any(axis=1)
# Drop rows with any missing data
penguins = penguins.loc[~missing_data]
penguins.shape

(333, 8)

Otteniamo in questo modo un DataFrame con 333 righe e 8 colonne. Calcoliamo la media del peso body_mass_g in funzione della specie:

summary_stats = (
    penguins.loc[:, ["species", "body_mass_g"]]
    .groupby("species")
    .aggregate(["mean", "std", "count"])
)
summary_stats.round(1)

	body_mass_g
	mean	std	count
species
Adelie	3706.2	458.6	146
Chinstrap	3733.1	384.3	68
Gentoo	5092.4	501.5	119

Creiamo un grafico a violino per questi dati.

_ = sns.catplot(kind="violin", data=penguins, x="species", y="body_mass_g")

Fino ad ora abbiamo analizzato un campione di osservazioni. Tuttavia, la stima della media e della varianza ottenuta dal campione non è rappresentativa dell’intera popolazione di pinguini. Ci chiediamo, quindi, quale sia l’incertezza associata alla stima del peso per le tre specie di pinguini, considerando che abbiamo a disposizione solo un piccolo campione di osservazioni.

Una possibile soluzione per quantificare l’incertezza della nostra stima è l’utilizzo dei metodi bayesiani.

Per fare ciò, dobbiamo definire un modello statistico che descriva la relazione tra i dati e i parametri. Nel nostro caso, descriviamo la distribuzione a posteriori dei parametri, tenendo conto dei dati, come segue:

\[ p(\mu, \sigma \mid y) \propto p(y \mid \mu, \sigma) p(\mu)p(\sigma), \]

dove la verosimiglianza sarà una densità Normale, dipendente da due parametri: \(\mu\) e \(\sigma\). Inoltre, per definire le distribuzioni a priori di questi due parametri, useremo una distribuzione Normale per \(\mu\) (\(\mu \sim \mathcal{N}(4000, 2000)\)) e una distribuzione Normale troncata per \(\sigma\) (valori negativi sono impossibili: \(\sigma \sim \mathcal{HN}(2000)\)).

Per esempio, selezioniamo solo i pinguini della specie Adelie. I valori della variabile body_mass_g relativi a questa specie sono:

adelie_mass_obs = penguins[penguins["species"] == "Adelie"]["body_mass_g"]
print(adelie_mass_obs)

0      3750.0
1      3800.0
2      3250.0
4      3450.0
5      3650.0
        ...  
147    3475.0
148    3450.0
149    3750.0
150    3700.0
151    4000.0
Name: body_mass_g, Length: 146, dtype: float64

Per implementare il modello descritto in precedenza, utilizzeremo PyMC.

Prima di eseguire il campionamento della distribuzione a posteriori, esamineremo la distribuzione predittiva a priori per verificare l’adeguatezza delle assunzioni relative alle distribuzioni a priori dei parametri. Per fare ciò, useremo la funzione pm.sample_prior_predictive(), che ci permetterà di generare 1000 valori casuali dalle distribuzioni a priori dei due parametri.

with pm.Model() as model_adelie_penguin_mass:
    sigma = pm.HalfNormal("sigma", sigma=1000)
    mu = pm.Normal("mu", 4000, 2000)
    mass = pm.Normal("mass", mu=mu, sigma=sigma, observed=adelie_mass_obs)

    idata1 = pm.sample_prior_predictive(samples=1000, random_seed=rng)

Sampling: [mass, mu, sigma]

Dall’oggetto inferenceData idata1 recuperiamo i campioni delle distribuzioni a priori e li passiamo a az.plot_trace():

prior = idata1.prior
_ = az.plot_trace(prior)
plt.tight_layout()

Dall’osservazione delle distribuzioni a priori dei due parametri, si può notare che non stiamo limitando troppo il campionamento della variabile adelie_mass_obs. Potremmo persino aver scelto distribuzioni a priori troppo ampie, dato che la distribuzione a priori per \(\mu\) include anche valori negativi. Tuttavia, poiché questo è un modello semplice e abbiamo un numero sufficiente di osservazioni, non c’è bisogno di preoccuparsi. Possiamo procedere con la stima della distribuzione a posteriori.

with pm.Model() as model_adelie_penguin_mass:
    
    sigma = pm.HalfNormal("sigma", sigma=1000)
    mu = pm.Normal("mu", 4000, 2000)
    mass = pm.Normal("mass", mu=mu, sigma=sigma, observed=adelie_mass_obs)

    idata2 = pm.sample(3000, chains=4)

Auto-assigning NUTS sampler...

Initializing NUTS using jitter+adapt_diag...

Multiprocess sampling (4 chains in 4 jobs)

NUTS: [sigma, mu]

---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
Cell In[9], line 7
      4 mu = pm.Normal("mu", 4000, 2000)
      5 mass = pm.Normal("mass", mu=mu, sigma=sigma, observed=adelie_mass_obs)
----> 7 idata2 = pm.sample(3000, chains=4)

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)
    764 _print_step_hierarchy(step)
    765 try:
--> 766     _mp_sample(**sample_args, **parallel_args)
    767 except pickle.PickleError:
    768     _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)
   1138 # We did draws += tune in pm.sample
   1139 draws -= tune
-> 1141 sampler = ps.ParallelSampler(
   1142     draws=draws,
   1143     tune=tune,
   1144     chains=chains,
   1145     cores=cores,
   1146     seeds=random_seed,
   1147     start_points=start,
   1148     step_method=step,
   1149     progressbar=progressbar,
   1150     mp_ctx=mp_ctx,
   1151 )
   1152 try:
   1153     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)
    399 if mp_ctx.get_start_method() != "fork":
    400     step_method_pickled = cloudpickle.dumps(step_method, protocol=-1)
--> 402 self._samplers = [
    403     ProcessAdapter(
    404         draws,
    405         tune,
    406         step_method,
    407         step_method_pickled,
    408         chain,
    409         seed,
    410         start,
    411         mp_ctx,
    412     )
    413     for chain, seed, start in zip(range(chains), seeds, start_points)
    414 ]
    416 self._inactive = self._samplers.copy()
    417 self._finished: List[ProcessAdapter] = []

File ~/mambaforge/envs/pymc_env/lib/python3.11/site-packages/pymc/sampling/parallel.py:403, in <listcomp>(.0)
    399 if mp_ctx.get_start_method() != "fork":
    400     step_method_pickled = cloudpickle.dumps(step_method, protocol=-1)
    402 self._samplers = [
--> 403     ProcessAdapter(
    404         draws,
    405         tune,
    406         step_method,
    407         step_method_pickled,
    408         chain,
    409         seed,
    410         start,
    411         mp_ctx,
    412     )
    413     for chain, seed, start in zip(range(chains), seeds, start_points)
    414 ]
    416 self._inactive = self._samplers.copy()
    417 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)
    242     step_method_send = step_method
    244 self._process = mp_ctx.Process(
    245     daemon=True,
    246     name=process_name,
   (...)
    257     ),
    258 )
--> 259 self._process.start()
    260 # Close the remote pipe, so that we get notified if the other
    261 # end is closed.
    262 remote_conn.close()

File ~/mambaforge/envs/pymc_env/lib/python3.11/multiprocessing/process.py:121, in BaseProcess.start(self)
    118 assert not _current_process._config.get('daemon'), \
    119        'daemonic processes are not allowed to have children'
    120 _cleanup()
--> 121 self._popen = self._Popen(self)
    122 self._sentinel = self._popen.sentinel
    123 # Avoid a refcycle if the target function holds an indirect
    124 # 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)
    297 @staticmethod
    298 def _Popen(process_obj):
    299     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)
     33 def __init__(self, process_obj):
     34     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)
     17 self.returncode = None
     18 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)
     55 self.finalizer = util.Finalize(self, util.close_fds,
     56                                (_parent_w, self.sentinel))
     57 with open(w, 'wb', closefd=True) as f:
---> 58     f.write(buf.getbuffer())
     59 self.pid = forkserver.read_signed(self.sentinel)

KeyboardInterrupt: 

pm.model_to_graphviz(model_adelie_penguin_mass)

Il KDE plot e il rank plot della distribuzione a posteriori dei parametri del modello bayesiano ci aiutano a diagnosticare eventuali problemi durante il campionamento in modo visivo. Tuttavia, in questo caso, non ci sono segnali di problemi durante il campionamento.

axes = az.plot_trace(idata2, divergences="bottom", kind="rank_bars")
plt.tight_layout()

Possiamo unire le quattro catene ottenute dal campionamento e utilizzarle per fare una stima più precisa dei parametri \(\mu\) e \(\sigma\) della distribuzione a posteriori.

Nella figura seguente possiamo vedere un confronto tra la distribuzione a posteriori stimata e le stime puntuali dei due parametri che abbiamo ottenuto dal campione. Le linee verticali indicano i valori puntuali stimati per \(\mu\) e \(\sigma\), mentre la curva rappresenta la distribuzione a posteriori stimata attraverso il campionamento delle catene.

axes = az.plot_posterior(idata2, hdi_prob=0.94)

axes[0].axvline(3706.2, linestyle="--")
axes[1].axvline(458.6, linestyle="--")

<matplotlib.lines.Line2D at 0x18c57fd50>

La stima bayesiana ci fornisce una distribuzione dei valori plausibili dei parametri. Le distribuzioni a posteriori ottenute per \(\mu\) e \(\sigma\) ci indicano il valore più credibile, ma anche la gamma di valori in cui ci si può aspettare che si trovi il “vero” valore del parametro, considerando differenti campioni di osservazioni. L’intervallo di valori plausibili è abbastanza stretto sia per \(\mu\) (tra 3634 e 3778), che per \(\sigma\) (94% HDI: [412, 513]).

In definitiva, le due distribuzioni a posteriori descrivono l’incertezza che abbiamo relativamente ai veri valori di \(\mu\) e \(\sigma\), tenendo in considerazione le nostre credenze precedenti (in questo caso, distribuzioni a priori debolmente informative, utilizzate per regolarizzare i dati), i dati osservati e l’ipotesi che abbiamo fatto sul meccanismo generatore dei dati.

Un riassunto numerico può essere ottenuto utilizzando le seguenti istruzioni:

az.summary(idata2).round(2)

	mean	sd	hdi_3%	hdi_97%	mcse_mean	mcse_sd	ess_bulk	ess_tail	r_hat
mu	3706.25	38.4	3634.82	3779.04	0.35	0.25	11858.0	8525.0	1.0
sigma	462.49	27.4	412.13	513.78	0.25	0.18	11766.0	8595.0	1.0

40.2. Confronto tra più gruppi

Per confrontare più gruppi, è necessario eseguire la stessa procedura di campionamento che abbiamo descritto precedentemente separatamente per ciascun gruppo di osservazioni. Pertanto, dobbiamo trovare un modo per comunicare a PyMC di effettuare questa operazione.

Prima di comprendere la procedura da seguire, dobbiamo capire come i dati sono organizzati in un oggetto della classe inferenceData.

40.3. Coordinate in PyMC e negli oggetti InferenceData

I risultati dell’analisi a posteriori ottenuti tramite PyMC vengono memorizzati in un oggetto di tipo inferenceData. A volte, per analizzare questi dati, è necessario estrarli dall’oggetto inferenceData. Per capire come ciò sia possibile, analizzeremo in dettaglio le proprietà dell’oggetto inferenceData.

L’oggetto inferenceData di ArviZ è in realtà un oggetto di tipo Dataset di xarray. xarray è una libreria che generalizza sia gli array di Numpy che i dataframe di Pandas. Di seguito viene mostrato un esempio di un array standard di Numpy.

npdata = rng.standard_normal(size=(2, 3))
npdata

array([[ 1.32565874, -0.57097002,  0.6682492 ],
       [ 1.14538102,  0.58734406, -0.05144329]])

Useremo questo array di Numpy per creare un oggetto DataArray di xarray.

data = xr.DataArray(
    npdata,
    dims=("user", "day"),
    coords={"user": ["Alice", "Bob"], "day": ["yesterday", "today", "tomorrow"]},
)
data

<xarray.DataArray (user: 2, day: 3)>
array([[ 1.32565874, -0.57097002,  0.6682492 ],
       [ 1.14538102,  0.58734406, -0.05144329]])
Coordinates:
  * user     (user) <U5 'Alice' 'Bob'
  * day      (day) <U9 'yesterday' 'today' 'tomorrow'

xarray.DataArray

• user: 2
• day: 3

• 1.326 -0.571 0.6682 1.145 0.5873 -0.05144

  array([[ 1.32565874, -0.57097002,  0.6682492 ],
         [ 1.14538102,  0.58734406, -0.05144329]])

• Coordinates: (2)
  • user
    (user)
    <U5
    'Alice' 'Bob'

    array(['Alice', 'Bob'], dtype='<U5')

  • day
    (day)
    <U9
    'yesterday' 'today' 'tomorrow'

    array(['yesterday', 'today', 'tomorrow'], dtype='<U9')

• Indexes: (2)
  • user
    PandasIndex

    PandasIndex(Index(['Alice', 'Bob'], dtype='object', name='user'))

  • day
    PandasIndex

    PandasIndex(Index(['yesterday', 'today', 'tomorrow'], dtype='object', name='day'))

• Attributes: (0)

Nell’esempio mostrato, i dati non sono stati modificati ma sono stati arricchiti con informazioni aggiuntive. L’oggetto xarray è di classe DataArray e rappresenta un array con due dimensioni: una chiamata “user” e l’altra chiamata “day”.

data.dims

('user', 'day')

Ogni dimensione è composta da diverse coordinate. La dimensione “user” comprende le coordinate “Alice” e “Bob”, mentre la dimensione “day” include le coordinate “yesterday”, “today” e “tomorrow”.

data.coords

Coordinates:
  * user     (user) <U5 'Alice' 'Bob'
  * day      (day) <U9 'yesterday' 'today' 'tomorrow'

Le coordinate in xarray possono essere utilizzate al posto degli indici numerici standard. Ad esempio, per selezionare il valore corrispondente a “Alice” e “yesterday” dall’array di Numpy, possiamo usare le coordinate nel seguente modo.

npdata[0,0]

1.3256587436158882

Eseguiamo nuovamente la stessa operazione usando questa volta l’oggetto DataArray di xarray.

data.sel(user="Alice", day="yesterday")

<xarray.DataArray ()>
array(1.32565874)
Coordinates:
    user     <U5 'Alice'
    day      <U9 'yesterday'

xarray.DataArray

• 1.326

  array(1.32565874)

• Coordinates: (2)
  • user
    ()
    <U5
    'Alice'

    array('Alice', dtype='<U5')

  • day
    ()
    <U9
    'yesterday'

    array('yesterday', dtype='<U9')

• Indexes: (0)
• Attributes: (0)

Uno dei vantaggi degli oggetti xarray consiste nella presenza di coordinate semanticamente etichettate. Non è più necessario tenere a mente che npdata[:,0] corrisponde a “yesterday” e npdata[1,:] corrisponde a “Bob”. Inoltre, non è nemmeno necessario ricordare l’ordine degli assi. L’ordine in cui le dimensioni sono memorizzate è irrilevante poiché usiamo le etichette delle dimensioni. Questo rende l’indicizzazione più intuitiva e meno soggetta a errori.

data.sel(day="yesterday", user="Alice")

<xarray.DataArray ()>
array(1.32565874)
Coordinates:
    user     <U5 'Alice'
    day      <U9 'yesterday'

xarray.DataArray

• 1.326

  array(1.32565874)

• Coordinates: (2)
  • user
    ()
    <U5
    'Alice'

    array('Alice', dtype='<U5')

  • day
    ()
    <U9
    'yesterday'

    array('yesterday', dtype='<U9')

• Indexes: (0)
• Attributes: (0)

data.sel(user="Alice", day="yesterday")

<xarray.DataArray ()>
array(1.32565874)
Coordinates:
    user     <U5 'Alice'
    day      <U9 'yesterday'

xarray.DataArray

• 1.326

  array(1.32565874)

• Coordinates: (2)
  • user
    ()
    <U5
    'Alice'

    array('Alice', dtype='<U5')

  • day
    ()
    <U9
    'yesterday'

    array('yesterday', dtype='<U9')

• Indexes: (0)
• Attributes: (0)

Inoltre, xarray permette la creazione di oggetti contenenti più oggetti DataArray, che prendono il nome di Dataset. Quando si ottiene un oggetto inferenceData di ArviZ dopo il campionamento, si ha a che fare con un oggetto Dataset di xarray. Consdieriamo un esempio.

var2 = xr.DataArray(
    rng.standard_normal(size=(2, 2)),
    dims=("x", "y"),
    coords={"x": [0, 1], "y": [11, 42]},
)
var3 = xr.DataArray(
    rng.standard_normal(size=(2, 2)),
    dims=("a", "b"),
    coords={"a": [4.2, 11.8], "b": ["Geneva", "London"]},
)
ds = xr.Dataset(dict(orig=data, v2=var2, v3=var3))
ds

<xarray.Dataset>
Dimensions:  (user: 2, day: 3, x: 2, y: 2, a: 2, b: 2)
Coordinates:
  * user     (user) <U5 'Alice' 'Bob'
  * day      (day) <U9 'yesterday' 'today' 'tomorrow'
  * x        (x) int64 0 1
  * y        (y) int64 11 42
  * a        (a) float64 4.2 11.8
  * b        (b) <U6 'Geneva' 'London'
Data variables:
    orig     (user, day) float64 1.326 -0.571 0.6682 1.145 0.5873 -0.05144
    v2       (x, y) float64 -1.845 -0.3072 -0.01176 0.6354
    v3       (a, b) float64 0.1133 1.419 0.7994 -0.4276

xarray.Dataset

• Dimensions:
  • user: 2
  • day: 3
  • x: 2
  • y: 2
  • a: 2
  • b: 2
• Coordinates: (6)
  • user
    (user)
    <U5
    'Alice' 'Bob'

    array(['Alice', 'Bob'], dtype='<U5')

  • day
    (day)
    <U9
    'yesterday' 'today' 'tomorrow'

    array(['yesterday', 'today', 'tomorrow'], dtype='<U9')

  • x
    (x)
    int64
    0 1

    array([0, 1])

  • y
    (y)
    int64
    11 42

    array([11, 42])

  • a
    (a)
    float64
    4.2 11.8

    array([ 4.2, 11.8])

  • b
    (b)
    <U6
    'Geneva' 'London'

    array(['Geneva', 'London'], dtype='<U6')

• Data variables: (3)
  • orig
    (user, day)
    float64
    1.326 -0.571 ... 0.5873 -0.05144

    array([[ 1.32565874, -0.57097002,  0.6682492 ],
           [ 1.14538102,  0.58734406, -0.05144329]])

  • v2
    (x, y)
    float64
    -1.845 -0.3072 -0.01176 0.6354

    array([[-1.84504794, -0.30722036],
           [-0.01176056,  0.63543926]])

  • v3
    (a, b)
    float64
    0.1133 1.419 0.7994 -0.4276

    array([[ 0.11328912,  1.41879627],
           [ 0.79941587, -0.42757956]])

• Indexes: (6)
  • user
    PandasIndex

    PandasIndex(Index(['Alice', 'Bob'], dtype='object', name='user'))

  • day
    PandasIndex

    PandasIndex(Index(['yesterday', 'today', 'tomorrow'], dtype='object', name='day'))

  • x
    PandasIndex

    PandasIndex(Int64Index([0, 1], dtype='int64', name='x'))

  • y
    PandasIndex

    PandasIndex(Int64Index([11, 42], dtype='int64', name='y'))

  • a
    PandasIndex

    PandasIndex(Float64Index([4.2, 11.8], dtype='float64', name='a'))

  • b
    PandasIndex

    PandasIndex(Index(['Geneva', 'London'], dtype='object', name='b'))

• Attributes: (0)

Si noti che il Dataset ha 6 dimensioni (user, day, x, y, a, b), ma ogni variabile all’interno del Dataset utilizza solo 2 di queste dimensioni.

Per accedere a ogni variabile singolarmente, possiamo utilizzare l’indicizzazione dei dizionari:

ds["orig"]

<xarray.DataArray 'orig' (user: 2, day: 3)>
array([[ 1.32565874, -0.57097002,  0.6682492 ],
       [ 1.14538102,  0.58734406, -0.05144329]])
Coordinates:
  * user     (user) <U5 'Alice' 'Bob'
  * day      (day) <U9 'yesterday' 'today' 'tomorrow'

xarray.DataArray

'orig'

• user: 2
• day: 3

• 1.326 -0.571 0.6682 1.145 0.5873 -0.05144

  array([[ 1.32565874, -0.57097002,  0.6682492 ],
         [ 1.14538102,  0.58734406, -0.05144329]])

• Coordinates: (2)
  • user
    (user)
    <U5
    'Alice' 'Bob'

    array(['Alice', 'Bob'], dtype='<U5')

  • day
    (day)
    <U9
    'yesterday' 'today' 'tomorrow'

    array(['yesterday', 'today', 'tomorrow'], dtype='<U9')

• Indexes: (2)
  • user
    PandasIndex

    PandasIndex(Index(['Alice', 'Bob'], dtype='object', name='user'))

  • day
    PandasIndex

    PandasIndex(Index(['yesterday', 'today', 'tomorrow'], dtype='object', name='day'))

• Attributes: (0)

ds["orig"].sel(user="Alice")

<xarray.DataArray 'orig' (day: 3)>
array([ 1.32565874, -0.57097002,  0.6682492 ])
Coordinates:
    user     <U5 'Alice'
  * day      (day) <U9 'yesterday' 'today' 'tomorrow'

xarray.DataArray

'orig'

• day: 3

• 1.326 -0.571 0.6682

  array([ 1.32565874, -0.57097002,  0.6682492 ])

• Coordinates: (2)
  • user
    ()
    <U5
    'Alice'

    array('Alice', dtype='<U5')

  • day
    (day)
    <U9
    'yesterday' 'today' 'tomorrow'

    array(['yesterday', 'today', 'tomorrow'], dtype='<U9')

• Indexes: (1)
  • day
    PandasIndex

    PandasIndex(Index(['yesterday', 'today', 'tomorrow'], dtype='object', name='day'))

• Attributes: (0)

ds["orig"].sel(user="Alice", day="yesterday")

<xarray.DataArray 'orig' ()>
array(1.32565874)
Coordinates:
    user     <U5 'Alice'
    day      <U9 'yesterday'

xarray.DataArray

'orig'

• 1.326

  array(1.32565874)

• Coordinates: (2)
  • user
    ()
    <U5
    'Alice'

    array('Alice', dtype='<U5')

  • day
    ()
    <U9
    'yesterday'

    array('yesterday', dtype='<U9')

• Indexes: (0)
• Attributes: (0)

40.4. Confonto tra gruppi

Analizziamo ora un esempio che ci permetterà di applicare le nozioni discusse in precedenza per effettuare il campionamento in modo separato su diversi gruppi di osservazioni, così da poter eseguire l’inferenza riguardante le differenze tra le medie dei gruppi.

Per indicare l’appartenenza al gruppo attraverso degli indici numerici, utilizziamo la funzione pd.factorize per creare dei valori numerici distinti (0, 1, 2) che corrispondono ai tre gruppi. Questo indice numerico verrà chiamato species_idx.

species_idx, species_codes = pd.factorize(penguins["species"])
print(species_idx)

[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1
 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2]

Il valore 0 identifica la specie Amelie, il valore 1 identifica la specie Gentoo e il valore 2 la specie Chinstrap.

print(*penguins["species"])

Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Gentoo Chinstrap Chinstrap Chinstrap Chinstrap Chinstrap Chinstrap Chinstrap Chinstrap Chinstrap Chinstrap Chinstrap Chinstrap Chinstrap Chinstrap Chinstrap Chinstrap Chinstrap Chinstrap Chinstrap Chinstrap Chinstrap Chinstrap Chinstrap Chinstrap Chinstrap Chinstrap Chinstrap Chinstrap Chinstrap Chinstrap Chinstrap Chinstrap Chinstrap Chinstrap Chinstrap Chinstrap Chinstrap Chinstrap Chinstrap Chinstrap Chinstrap Chinstrap Chinstrap Chinstrap Chinstrap Chinstrap Chinstrap Chinstrap Chinstrap Chinstrap Chinstrap Chinstrap Chinstrap Chinstrap Chinstrap Chinstrap Chinstrap Chinstrap Chinstrap Chinstrap Chinstrap Chinstrap Chinstrap Chinstrap Chinstrap Chinstrap Chinstrap Chinstrap

Di seguito è riportata la descrizione del modello per il campionamento MCMC che distingue tra più gruppi di osservazioni:

Per prima cosa, utilizziamo pd.factorize per creare dei valori numerici che distinguono tra i tre gruppi di osservazioni, chiamati species_idx. Successivamente, specifichiamo la distribuzione a priori per mu e sigma per ciascun gruppo di osservazioni. In questo caso, abbiamo aggiunto l’argomento dims = "species" nella specificazione delle distribuzioni a priori. Nella verosimiglianza, usiamo mu_prior[species_idx] e sigma_prior[species_idx] per indicare che le osservazioni appartengono a tre gruppi diversi, come specificato dalla variabile species. Inoltre, passiamo il dizionario coords a pm.Model(), dove la chiave species è associata ad una lista di interi [0, 1, 2]. Il modello per il campionamento MCMC che distingue tra più gruppi di osservazioni ha quindi una specifica quasi identica a quella usata per un singolo gruppo di osservazioni, con la sola aggiunta della specificazione dei dims e dei valori numerici associati ai gruppi di osservazioni.

coords = {"species": [0, 1, 2]}

with pm.Model(coords=coords) as labeled_model:
    mu_prior = pm.Normal("mu", mu=4000, sigma=500, dims="species")
    sigma_prior = pm.HalfCauchy("sigma", beta=1000, dims="species")
    likelihood = pm.Normal(
        "likelihood",
        mu=mu_prior[species_idx],
        sigma=sigma_prior[species_idx],
        observed=penguins["body_mass_g"],
    )

Eseguiamo il campionamento.

with labeled_model:
    labeled_idata = pm.sample()

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

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

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

Esaminiamo le distribuzioni a posteriori e le diagnostiche del campionamento.

axes = az.plot_trace(
    labeled_idata, combined=True, divergences="bottom", kind="rank_bars"
)
plt.tight_layout()

Abbiamo ottenuto tre distribuzioni a posteriori per ciascuno dei due parametri \(\mu\) e \(\sigma\), una distribuzione per ciascuna specie di pinguini.

Esaminiamo se le quattro catene hanno prodotto risultati diversi per il parametro \(\mu\).

axes = az.plot_forest(labeled_idata, var_names=["mu"], figsize=(8, 3))
axes[0].set_title("μ Mass Estimate: 94.0% HDI")

Text(0.5, 1.0, 'μ Mass Estimate: 94.0% HDI')

Facciamo la stessa cosa per il parametro \(\sigma\).

axes = az.plot_forest(labeled_idata, var_names=["sigma"], figsize=(8, 3))
axes[0].set_title("σ Mass Estimate: 94.0% HDI")

Text(0.5, 1.0, 'σ Mass Estimate: 94.0% HDI')

Un sommario numerico delle distribuzioni a posteriori si ottiene nel modo seguente.

az.summary(labeled_idata, kind="stats").round(2)

	mean	sd	hdi_3%	hdi_97%
mu[0]	3707.94	37.59	3641.64	3779.60
mu[1]	5083.49	45.63	4994.39	5164.83
mu[2]	3736.14	46.50	3650.66	3826.96
sigma[0]	461.56	27.02	408.23	511.25
sigma[1]	505.64	33.34	448.57	572.55
sigma[2]	390.13	34.40	328.97	455.25

Dalle distribuzioni a posteriori e dagli intervalli di credibilità, possiamo dedurre che il peso dei pinguini Adelie e Chinstrap è simile, mentre i pinguini Gentoo sono in genere più pesanti. Inoltre, l’analisi indica che il peso dei pinguini Gentoo ha una maggiore variabilità rispetto ai pinguini Adelie, seguita da quella dei pinguini Chinstrap. Al contrario, i pinguini Chinstrap mostrano una maggiore omogeneità nel peso.

40.5. Verifica di ipotesi bayesiana

Una volta ottenuto un campione della distribuzione a priori del parametro di interesse \(\mu\) per ciascuna delle tre specie di pinguini, possiamo chiederci quale sia la probabilità che il peso di un pinguino di una specie sia maggiore di quello di un pinguino di un’altra specie. Possiamo rispondere a questa domanda utilizzando un campione casuale della distribuzione a posteriori. Per fare ciò, dobbiamo confrontare il parametro di interesse per molti valori e calcolare la media dei valori trovati.

Per recuperare i valori numerici della distribuzione a posteriori del parametro \(\mu\) dall’oggetto labeled_idata, possiamo utilizzare le funzioni di ArviZ. Esaminiamo le proprietà dell’oggetto labeled_idata.

labeled_idata

arviz.InferenceData

• posterior

  <xarray.Dataset>
  Dimensions:  (chain: 4, draw: 1000, species: 3)
  Coordinates:
    * chain    (chain) int64 0 1 2 3
    * draw     (draw) int64 0 1 2 3 4 5 6 7 8 ... 992 993 994 995 996 997 998 999
    * species  (species) int64 0 1 2
  Data variables:
      mu       (chain, draw, species) float64 3.713e+03 5.1e+03 ... 3.738e+03
      sigma    (chain, draw, species) float64 442.9 475.1 381.8 ... 515.7 370.2
  Attributes:
      created_at:                 2023-05-06T06:30:02.078800
      arviz_version:              0.15.1
      inference_library:          pymc
      inference_library_version:  5.3.0
      sampling_time:              43.66927099227905
      tuning_steps:               1000

  xarray.Dataset

  • Dimensions:
    • chain: 4
    • draw: 1000
    • species: 3
  • Coordinates: (3)
    • chain
      (chain)
      int64
      0 1 2 3

      array([0, 1, 2, 3])

    • draw
      (draw)
      int64
      0 1 2 3 4 5 ... 995 996 997 998 999

      array([  0,   1,   2, ..., 997, 998, 999])

    • species
      (species)
      int64
      0 1 2

      array([0, 1, 2])

  • Data variables: (2)
    • mu
      (chain, draw, species)
      float64
      3.713e+03 5.1e+03 ... 3.738e+03

      array([[[3713.24652133, 5099.87615717, 3798.15788846],
              [3694.75257683, 5085.88745042, 3738.96050014],
              [3721.56984315, 5077.98993634, 3707.89450852],
              ...,
              [3717.41427976, 5127.70598267, 3582.2105841 ],
              [3724.37072298, 5111.63450185, 3747.81900595],
              [3707.57508782, 5012.79650676, 3702.71402215]],
      
             [[3710.57791675, 5065.31239168, 3707.38523188],
              [3725.5608288 , 5122.7599556 , 3748.50483432],
              [3711.68885019, 5079.56815309, 3730.05125194],
              ...,
              [3713.41049886, 5137.02745921, 3598.34502641],
              [3638.51049066, 5130.76412867, 3692.56016114],
              [3774.78193344, 5061.6994933 , 3742.52586122]],
      
             [[3710.60894828, 5062.8548782 , 3645.76084736],
              [3651.85042772, 5111.99483128, 3665.79353687],
              [3700.60974818, 5066.96782879, 3641.39431661],
              ...,
              [3724.92663273, 5098.41457784, 3779.33694367],
              [3685.5930511 , 5038.16586981, 3743.17490689],
              [3731.25053034, 5106.53654003, 3680.40519101]],
      
             [[3660.58320382, 5087.97917745, 3716.43412258],
              [3708.67099564, 5090.15320339, 3763.62840167],
              [3709.23995653, 5078.73833181, 3715.92500142],
              ...,
              [3703.11953167, 5037.25969865, 3745.22845812],
              [3677.56865981, 5099.46450144, 3680.75493962],
              [3741.3153177 , 5062.55402125, 3738.07389079]]])

    • sigma
      (chain, draw, species)
      float64
      442.9 475.1 381.8 ... 515.7 370.2

      array([[[442.91054142, 475.11410774, 381.78278037],
              [496.20557262, 536.81825462, 404.41421935],
              [421.81269784, 475.7924003 , 366.7388128 ],
              ...,
              [462.26982112, 502.72798581, 364.30805191],
              [503.68637691, 491.32010969, 341.87734819],
              [497.29922032, 498.53840384, 391.15587744]],
      
             [[436.33876257, 464.28851134, 379.90316726],
              [493.61101205, 541.4704704 , 383.03562888],
              [464.77518935, 466.5708577 , 373.54371054],
              ...,
              [491.91439533, 462.38377871, 375.81745359],
              [471.3987055 , 537.79860339, 440.12169636],
              [446.29835339, 468.53731585, 363.034989  ]],
      
             [[498.69306497, 564.40319463, 515.32724557],
              [455.0173386 , 485.84666566, 345.15085577],
              [399.05449896, 550.21215501, 349.83453141],
              ...,
              [478.08336024, 482.3774175 , 357.12874992],
              [475.71698743, 491.87775996, 409.52843093],
              [505.46473166, 479.22302303, 425.40553756]],
      
             [[460.41954009, 467.87356662, 468.22100133],
              [407.18557044, 489.20046564, 447.38664368],
              [519.55053618, 501.44583546, 343.50864277],
              ...,
              [464.73415104, 505.24157962, 394.0705571 ],
              [477.28878831, 484.4172539 , 435.83634427],
              [441.43757043, 515.6576931 , 370.16304851]]])

  • Indexes: (3)
    • chain
      PandasIndex

      PandasIndex(Int64Index([0, 1, 2, 3], dtype='int64', name='chain'))

    • draw
      PandasIndex

      PandasIndex(Int64Index([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,
                  ...
                  990, 991, 992, 993, 994, 995, 996, 997, 998, 999],
                 dtype='int64', name='draw', length=1000))

    • species
      PandasIndex

      PandasIndex(Int64Index([0, 1, 2], dtype='int64', name='species'))

  • Attributes: (6)

    created_at :

    2023-05-06T06:30:02.078800

    arviz_version :

    0.15.1

    inference_library :

    pymc

    inference_library_version :

    5.3.0

    sampling_time :

    43.66927099227905

    tuning_steps :

    1000

  
  
• sample_stats

  <xarray.Dataset>
  Dimensions:                (chain: 4, draw: 1000)
  Coordinates:
    * chain                  (chain) int64 0 1 2 3
    * draw                   (draw) int64 0 1 2 3 4 5 ... 994 995 996 997 998 999
  Data variables: (12/17)
      diverging              (chain, draw) bool False False False ... False False
      lp                     (chain, draw) float64 -2.54e+03 ... -2.54e+03
      n_steps                (chain, draw) float64 3.0 7.0 3.0 3.0 ... 3.0 3.0 3.0
      tree_depth             (chain, draw) int64 2 3 2 2 2 3 2 2 ... 2 2 2 3 2 2 2
      max_energy_error       (chain, draw) float64 0.232 0.3705 ... 0.7087 1.87
      index_in_trajectory    (chain, draw) int64 -1 2 3 2 2 -1 ... 2 3 2 2 -1 -3
      ...                     ...
      reached_max_treedepth  (chain, draw) bool False False False ... False False
      step_size              (chain, draw) float64 0.9499 0.9499 ... 1.125 1.125
      energy                 (chain, draw) float64 2.541e+03 ... 2.546e+03
      perf_counter_diff      (chain, draw) float64 0.001301 0.002175 ... 0.000774
      smallest_eigval        (chain, draw) float64 nan nan nan nan ... nan nan nan
      energy_error           (chain, draw) float64 0.232 -0.01623 ... -0.4988
  Attributes:
      created_at:                 2023-05-06T06:30:02.095549
      arviz_version:              0.15.1
      inference_library:          pymc
      inference_library_version:  5.3.0
      sampling_time:              43.66927099227905
      tuning_steps:               1000

  xarray.Dataset

  • Dimensions:
    • chain: 4
    • draw: 1000
  • Coordinates: (2)
    • chain
      (chain)
      int64
      0 1 2 3

      array([0, 1, 2, 3])

    • draw
      (draw)
      int64
      0 1 2 3 4 5 ... 995 996 997 998 999

      array([  0,   1,   2, ..., 997, 998, 999])

  • Data variables: (17)
    • diverging
      (chain, draw)
      bool
      False False False ... False False

      array([[False, False, False, ..., False, False, False],
             [False, False, False, ..., False, False, False],
             [False, False, False, ..., False, False, False],
             [False, False, False, ..., False, False, False]])

    • lp
      (chain, draw)
      float64
      -2.54e+03 -2.54e+03 ... -2.54e+03

      array([[-2540.24567864, -2540.40638374, -2540.52072521, ...,
              -2545.47281784, -2541.29551209, -2541.1350665 ],
             [-2540.09849127, -2540.67343707, -2539.45019638, ...,
              -2545.59073155, -2543.01932588, -2541.40359076],
             [-2547.33132134, -2542.34437209, -2545.89175023, ...,
              -2540.1948987 , -2539.93616396, -2541.83502   ],
             [-2542.50765489, -2542.57081279, -2541.85675853, ...,
              -2539.34373777, -2541.03011378, -2539.64469174]])

    • n_steps
      (chain, draw)
      float64
      3.0 7.0 3.0 3.0 ... 7.0 3.0 3.0 3.0

      array([[3., 7., 3., ..., 3., 7., 3.],
             [3., 7., 3., ..., 3., 3., 3.],
             [7., 7., 3., ..., 3., 3., 3.],
             [3., 7., 3., ..., 3., 3., 3.]])

    • tree_depth
      (chain, draw)
      int64
      2 3 2 2 2 3 2 2 ... 2 2 2 2 3 2 2 2

      array([[2, 3, 2, ..., 2, 3, 2],
             [2, 3, 2, ..., 2, 2, 2],
             [3, 3, 2, ..., 2, 2, 2],
             [2, 3, 2, ..., 2, 2, 2]])

    • max_energy_error
      (chain, draw)
      float64
      0.232 0.3705 0.1256 ... 0.7087 1.87

      array([[ 0.2319789 ,  0.37050795,  0.12559172, ...,  1.07286644,
              -2.26115951,  0.31222501],
             [-0.19785321,  0.73437493, -0.25236447, ...,  0.97568364,
              -0.68021195, -0.45126754],
             [ 1.05478524, -0.89265872,  1.04005716, ..., -0.31926886,
               0.5435203 ,  0.40946243],
             [ 1.26678014,  0.18459724, -0.23789352, ..., -0.34158518,
               0.70868795,  1.86996417]])

    • index_in_trajectory
      (chain, draw)
      int64
      -1 2 3 2 2 -1 -3 ... 2 3 2 2 -1 -3

      array([[-1,  2,  3, ...,  2, -1, -2],
             [ 2, -3,  2, ..., -2, -2,  3],
             [-5, -2, -1, ..., -2,  1,  2],
             [ 3,  1,  3, ...,  2, -1, -3]])

    • acceptance_rate
      (chain, draw)
      float64
      0.8696 0.9171 ... 0.5826 0.718

      array([[0.86955589, 0.917139  , 0.96065829, ..., 0.78067557, 0.97000679,
              0.89976886],
             [0.99087607, 0.65668364, 1.        , ..., 0.61046315, 1.        ,
              0.9619897 ],
             [0.72832311, 1.        , 0.4259881 , ..., 1.        , 0.8182979 ,
              0.78503514],
             [0.76057911, 0.92509348, 0.96302562, ..., 0.96285373, 0.58255362,
              0.71804306]])

    • step_size_bar
      (chain, draw)
      float64
      0.9484 0.9484 ... 0.934 0.934

      array([[0.94842329, 0.94842329, 0.94842329, ..., 0.94842329, 0.94842329,
              0.94842329],
             [0.92804114, 0.92804114, 0.92804114, ..., 0.92804114, 0.92804114,
              0.92804114],
             [0.86670509, 0.86670509, 0.86670509, ..., 0.86670509, 0.86670509,
              0.86670509],
             [0.93396366, 0.93396366, 0.93396366, ..., 0.93396366, 0.93396366,
              0.93396366]])

    • process_time_diff
      (chain, draw)
      float64
      0.001304 0.002176 ... 0.000775

      array([[0.001304, 0.002176, 0.000717, ..., 0.000677, 0.001301, 0.000669],
             [0.00075 , 0.001917, 0.001499, ..., 0.000882, 0.000753, 0.000845],
             [0.002704, 0.001466, 0.000727, ..., 0.001608, 0.001533, 0.001563],
             [0.001735, 0.002839, 0.001613, ..., 0.000677, 0.000665, 0.000775]])

    • perf_counter_start
      (chain, draw)
      float64
      4.32e+04 4.32e+04 ... 4.32e+04

      array([[43199.08027814, 43199.08189963, 43199.08427312, ...,
              43200.66852661, 43200.66935343, 43200.67080899],
             [43198.90462087, 43198.90555787, 43198.90771202, ...,
              43200.61216131, 43200.61330579, 43200.61452802],
             [43198.36011248, 43198.36334852, 43198.36501   , ...,
              43200.17638037, 43200.17953168, 43200.18165346],
             [43200.24358119, 43200.24586927, 43200.24911842, ...,
              43201.50485643, 43201.50569006, 43201.50651334]])

    • largest_eigval
      (chain, draw)
      float64
      nan nan nan nan ... nan nan nan nan

      array([[nan, nan, nan, ..., nan, nan, nan],
             [nan, nan, nan, ..., nan, nan, nan],
             [nan, nan, nan, ..., nan, nan, nan],
             [nan, nan, nan, ..., nan, nan, nan]])

    • reached_max_treedepth
      (chain, draw)
      bool
      False False False ... False False

      array([[False, False, False, ..., False, False, False],
             [False, False, False, ..., False, False, False],
             [False, False, False, ..., False, False, False],
             [False, False, False, ..., False, False, False]])

    • step_size
      (chain, draw)
      float64
      0.9499 0.9499 ... 1.125 1.125

      array([[0.94994045, 0.94994045, 0.94994045, ..., 0.94994045, 0.94994045,
              0.94994045],
             [1.07169211, 1.07169211, 1.07169211, ..., 1.07169211, 1.07169211,
              1.07169211],
             [0.8327984 , 0.8327984 , 0.8327984 , ..., 0.8327984 , 0.8327984 ,
              0.8327984 ],
             [1.12456665, 1.12456665, 1.12456665, ..., 1.12456665, 1.12456665,
              1.12456665]])

    • energy
      (chain, draw)
      float64
      2.541e+03 2.542e+03 ... 2.546e+03

      array([[2540.84668651, 2541.75664436, 2541.79372455, ..., 2546.96322596,
              2546.66173335, 2543.69950807],
             [2542.19895036, 2544.37734746, 2540.7542396 , ..., 2547.20286411,
              2548.24925405, 2545.4769314 ],
             [2553.20187791, 2548.98888028, 2550.15072615, ..., 2542.16275003,
              2542.59532279, 2543.09619631],
             [2548.11587879, 2545.70421969, 2545.1090972 , ..., 2541.64573858,
              2542.81826142, 2545.94327878]])

    • perf_counter_diff
      (chain, draw)
      float64
      0.001301 0.002175 ... 0.000774

      array([[0.00130128, 0.00217514, 0.0007166 , ..., 0.00067609, 0.00129972,
              0.00066903],
             [0.00074934, 0.00192281, 0.00155961, ..., 0.00092394, 0.00075111,
              0.00084565],
             [0.00300193, 0.00146488, 0.00072628, ..., 0.00163728, 0.00153227,
              0.00160638],
             [0.00195721, 0.00289511, 0.00202122, ..., 0.00067691, 0.00066485,
              0.00077401]])

    • smallest_eigval
      (chain, draw)
      float64
      nan nan nan nan ... nan nan nan nan

      array([[nan, nan, nan, ..., nan, nan, nan],
             [nan, nan, nan, ..., nan, nan, nan],
             [nan, nan, nan, ..., nan, nan, nan],
             [nan, nan, nan, ..., nan, nan, nan]])

    • energy_error
      (chain, draw)
      float64
      0.232 -0.01623 ... 0.3505 -0.4988

      array([[ 0.2319789 , -0.01622652, -0.00257769, ...,  1.07286644,
              -2.26115951,  0.03305046],
             [-0.19785321,  0.22986452, -0.25236447, ...,  0.97568364,
              -0.55766353, -0.45126754],
             [ 0.39060088, -0.82089416,  0.81964866, ..., -0.23699855,
               0.05886842,  0.3164157 ],
             [-0.11903678, -0.01554294, -0.23789352, ..., -0.34158518,
               0.35053153, -0.49883526]])

  • Indexes: (2)
    • chain
      PandasIndex

      PandasIndex(Int64Index([0, 1, 2, 3], dtype='int64', name='chain'))

    • draw
      PandasIndex

      PandasIndex(Int64Index([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,
                  ...
                  990, 991, 992, 993, 994, 995, 996, 997, 998, 999],
                 dtype='int64', name='draw', length=1000))

  • Attributes: (6)

    created_at :

    2023-05-06T06:30:02.095549

    arviz_version :

    0.15.1

    inference_library :

    pymc

    inference_library_version :

    5.3.0

    sampling_time :

    43.66927099227905

    tuning_steps :

    1000

  
  
• observed_data

  <xarray.Dataset>
  Dimensions:           (likelihood_dim_0: 333)
  Coordinates:
    * likelihood_dim_0  (likelihood_dim_0) int64 0 1 2 3 4 ... 328 329 330 331 332
  Data variables:
      likelihood        (likelihood_dim_0) float64 3.75e+03 3.8e+03 ... 3.775e+03
  Attributes:
      created_at:                 2023-05-06T06:30:02.102522
      arviz_version:              0.15.1
      inference_library:          pymc
      inference_library_version:  5.3.0

  xarray.Dataset

  • Dimensions:
    • likelihood_dim_0: 333
  • Coordinates: (1)
    • likelihood_dim_0
      (likelihood_dim_0)
      int64
      0 1 2 3 4 5 ... 328 329 330 331 332

      array([  0,   1,   2, ..., 330, 331, 332])

  • Data variables: (1)
    • likelihood
      (likelihood_dim_0)
      float64
      3.75e+03 3.8e+03 ... 3.775e+03

      array([3750., 3800., 3250., 3450., 3650., 3625., 4675., 3200., 3800.,
             4400., 3700., 3450., 4500., 3325., 4200., 3400., 3600., 3800.,
             3950., 3800., 3800., 3550., 3200., 3150., 3950., 3250., 3900.,
             3300., 3900., 3325., 4150., 3950., 3550., 3300., 4650., 3150.,
             3900., 3100., 4400., 3000., 4600., 3425., 3450., 4150., 3500.,
             4300., 3450., 4050., 2900., 3700., 3550., 3800., 2850., 3750.,
             3150., 4400., 3600., 4050., 2850., 3950., 3350., 4100., 3050.,
             4450., 3600., 3900., 3550., 4150., 3700., 4250., 3700., 3900.,
             3550., 4000., 3200., 4700., 3800., 4200., 3350., 3550., 3800.,
             3500., 3950., 3600., 3550., 4300., 3400., 4450., 3300., 4300.,
             3700., 4350., 2900., 4100., 3725., 4725., 3075., 4250., 2925.,
             3550., 3750., 3900., 3175., 4775., 3825., 4600., 3200., 4275.,
             3900., 4075., 2900., 3775., 3350., 3325., 3150., 3500., 3450.,
             3875., 3050., 4000., 3275., 4300., 3050., 4000., 3325., 3500.,
             3500., 4475., 3425., 3900., 3175., 3975., 3400., 4250., 3400.,
             3475., 3050., 3725., 3000., 3650., 4250., 3475., 3450., 3750.,
             3700., 4000., 4500., 5700., 4450., 5700., 5400., 4550., 4800.,
             5200., 4400., 5150., 4650., 5550., 4650., 5850., 4200., 5850.,
             4150., 6300., 4800., 5350., 5700., 5000., 4400., 5050., 5000.,
             5100., 5650., 4600., 5550., 5250., 4700., 5050., 6050., 5150.,
             5400., 4950., 5250., 4350., 5350., 3950., 5700., 4300., 4750.,
             5550., 4900., 4200., 5400., 5100., 5300., 4850., 5300., 4400.,
             5000., 4900., 5050., 4300., 5000., 4450., 5550., 4200., 5300.,
             4400., 5650., 4700., 5700., 5800., 4700., 5550., 4750., 5000.,
             5100., 5200., 4700., 5800., 4600., 6000., 4750., 5950., 4625.,
             5450., 4725., 5350., 4750., 5600., 4600., 5300., 4875., 5550.,
             4950., 5400., 4750., 5650., 4850., 5200., 4925., 4875., 4625.,
             5250., 4850., 5600., 4975., 5500., 5500., 4700., 5500., 4575.,
             5500., 5000., 5950., 4650., 5500., 4375., 5850., 6000., 4925.,
             4850., 5750., 5200., 5400., 3500., 3900., 3650., 3525., 3725.,
             3950., 3250., 3750., 4150., 3700., 3800., 3775., 3700., 4050.,
             3575., 4050., 3300., 3700., 3450., 4400., 3600., 3400., 2900.,
             3800., 3300., 4150., 3400., 3800., 3700., 4550., 3200., 4300.,
             3350., 4100., 3600., 3900., 3850., 4800., 2700., 4500., 3950.,
             3650., 3550., 3500., 3675., 4450., 3400., 4300., 3250., 3675.,
             3325., 3950., 3600., 4050., 3350., 3450., 3250., 4050., 3800.,
             3525., 3950., 3650., 3650., 4000., 3400., 3775., 4100., 3775.])

  • Indexes: (1)
    • likelihood_dim_0
      PandasIndex

      PandasIndex(Int64Index([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,
                  ...
                  323, 324, 325, 326, 327, 328, 329, 330, 331, 332],
                 dtype='int64', name='likelihood_dim_0', length=333))

  • Attributes: (4)

    created_at :

    2023-05-06T06:30:02.102522

    arviz_version :

    0.15.1

    inference_library :

    pymc

    inference_library_version :

    5.3.0

  
  

Recuperiamo i valori a posteriori del parametro \(\mu\).

mu_post = labeled_idata.posterior['mu']

L’oggett ottenuto è un array di dimensioni \(4 \times 1000 \times 3\). L’indice 4 si riferisce alle catente, l’indice 1000 si riferisce al numero di campioni ottenuti, l’indice 3 si riferisce ai gruppi.

mu_post.shape

(4, 1000, 3)

Per trovare la media a posteriori del parametro \(\mu\) della specie Adelie prendiamo le osservazioni che si riferiscono a tutte le catene (:), a tutti i valori campionati (:) per il primo gruppo (0). La media a posteriori dei pinguini Adelie è dunque uguale a

print(mu_post[:, :, 0].mean())

<xarray.DataArray 'mu' ()>
array(3707.93884115)
Coordinates:
    species  int64 0

Per tutte e tre le specie (Adelie, Gentoo, Chinstrap) abbiamo

np.array(
    [mu_post[:, :, 0].mean(), mu_post[:, :, 1].mean(), mu_post[:, :, 2].mean()]
).round(1)

array([3707.9, 5083.5, 3736.1])

summary_stats = (
    penguins.loc[:, ["species", "body_mass_g"]]
    .groupby("species")
    .aggregate(["mean", "std", "count"])
)
summary_stats.round(1)

	body_mass_g
	mean	std	count
species
Adelie	3706.2	458.6	146
Chinstrap	3733.1	384.3	68
Gentoo	5092.4	501.5	119

Per verificare che l’ordinamento dei gruppi corrisponda a Adelie, Gentoo, Chinstrap, esaminiamo la stima a posteriori della deviazione standard delle tre specie.

sigma_post = labeled_idata.posterior["sigma"]
np.array(
    [sigma_post[:, :, 0].mean(), sigma_post[:, :, 1].mean(), sigma_post[:, :, 2].mean()]
).round(1)

array([461.6, 505.6, 390.1])

I valori trovati confermano quello che ci aspettavamo. Ora che abbiamo capito come estrarre le stime a posteriori di un parametro per ciascun gruppo, possiamo passare al test di ipotesi bayesiano. Chiediamoci quale sia la probabilità che un pinguino Adelie abbia un peso minore di un pinguino Chinstrap.

np.mean([mu_post[:, :, 0] < mu_post[:, :, 2]])

0.68125

Chiediamoci quale sia la probabilità che un pinguino Adelie abbia un peso minore di un pinguino Gentoo.

np.mean([mu_post[:, :, 0] < mu_post[:, :, 2]])

0.68125

Chiediamoci quale sia la probabilità che un pinguino Adelie abbia un peso minore di un pinguino Gentoo.

np.mean([mu_post[:, :, 0] < mu_post[:, :, 1]])

1.0

Consideriamo ora le deviazioni standard. Troviamo la probabilità che la varianza del peso dei pinguini Adelie sia minore di quella dei pinguini Chinstrap.

np.mean([sigma_post[:, :, 2] > sigma_post[:, :, 0]])

0.055

40.6. Watermark

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

Last updated: Sat May 06 2023

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

pytensor: 2.11.1

pandas    : 2.0.1
pymc      : 5.3.0
numpy     : 1.24.3
matplotlib: 3.7.1
seaborn   : 0.12.2
xarray    : 2023.4.2
arviz     : 0.15.1

Watermark: 2.3.1