Metodi di sintesi della distribuzione a posteriori

L’obiettivo di questo capitolo è analizzare in dettaglio il flusso di lavoro all’interno del contesto dell’inferenza bayesiana, una volta che sia stata ottenuta la stima della distribuzione a posteriori mediante l’utilizzo della tecnica MCMC.

Preparazione del Notebook

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import pymc as pm
import pymc.sampling_jax
import arviz as az
import scipy.stats as stats
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")

Sintesi della Distribuzione a Posteriori

Uno dei fattori principali che ha contribuito alla crescente popolarità dei metodi bayesiani nelle scienze sociali è stata la (ri)scoperta di algoritmi numerici capaci di stimare le distribuzioni a posteriori dei parametri del modello a partire dai dati osservati. Prima di tali sviluppi, era praticamente impossibile ottenere misure riassuntive delle distribuzioni a posteriori, specialmente per modelli complessi con un elevato numero di parametri.

Gli algoritmi numerici che esamineremo in questo capitolo fanno uso dell’integrazione Monte Carlo attraverso catene di Markov, meglio conosciuta come campionamento Markov chain Monte Carlo (MCMC). Questa tecnica ha rivoluzionato la capacità di effettuare inferenze bayesiane, rendendo più accessibile e gestibile l’analisi di modelli complessi.

Campionamento con PyMC

A titolo di esempio, useremo ancora una volta il set di dati “moma_sample”, che è stato analizzato nel capitolo precedente. Iniziamo replicando il processo di campionamento PyMC descritto in precedenza.

y = 14
ntrials = 100

alpha_prior = 4
beta_prior = 6

with pm.Model() as bb_model:
    theta = pm.Beta("theta", alpha=alpha_prior, beta=beta_prior)
    obs = pm.Binomial("obs", p=theta, n=ntrials, observed=y)

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

Show code cell output Hide code cell output

Compiling...

Compilation time = 0:00:00.531103

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 0:   0%|                                                                 | 0/2000 [00:01<?, ?it/s]

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

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

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

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

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

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

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

Sampling time = 0:00:01.293610

Transforming variables...

Transformation time = 0:00:00.053729

La figura seguente rappresenta la distribuzione a posteriori di \(\theta\) ottenuta mediante il metodo MCMC e la traccia delle catene di Markov.

az.plot_trace(idata, combined=True)
plt.show()

Estrazione degli attributi da InferenceData

Ora che abbiamo ottenuto la distribuzione a posteriori, possiamo effettuare varie operazioni con essa:

• Stima puntuale

• Intervalli di credibilità

• Test di ipotesi

• Predizione

Stima puntuale

Procediamo all’esame dell’oggetto idata.

idata

arviz.InferenceData

• posterior

  <xarray.Dataset>
  Dimensions:  (chain: 4, draw: 1000)
  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
  Data variables:
      theta    (chain, draw) float64 0.1733 0.1732 0.2063 ... 0.1187 0.1444 0.1424
  Attributes:
      created_at:     2024-01-23T07:45:24.000626
      arviz_version:  0.17.0

  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: (1)
    • theta
      (chain, draw)
      float64
      0.1733 0.1732 ... 0.1444 0.1424

      array([[0.17326516, 0.17320269, 0.20629908, ..., 0.12454636, 0.16953616,
              0.18556175],
             [0.20671803, 0.14461032, 0.17306915, ..., 0.17527697, 0.1427207 ,
              0.15351922],
             [0.14086973, 0.19343519, 0.19594644, ..., 0.1142225 , 0.17836103,
              0.14123355],
             [0.11920471, 0.14753389, 0.26728206, ..., 0.11865411, 0.14442225,
              0.14238198]])

  • Indexes: (2)
    • chain
      PandasIndex

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

    • draw
      PandasIndex

      PandasIndex(Index([  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: (2)

    created_at :

    2024-01-23T07:45:24.000626

    arviz_version :

    0.17.0

  
  
• 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 6 ... 993 994 995 996 997 998 999
  Data variables:
      acceptance_rate  (chain, draw) float64 0.9147 0.7716 ... 0.6239 0.9669
      step_size        (chain, draw) float64 1.183 1.183 1.183 ... 0.9726 0.9726
      diverging        (chain, draw) bool False False False ... False False False
      energy           (chain, draw) float64 5.58 5.612 5.119 ... 5.44 9.76 4.973
      n_steps          (chain, draw) int64 3 3 3 1 3 1 3 1 3 ... 3 3 3 3 3 3 1 3 3
      tree_depth       (chain, draw) int64 2 2 2 1 2 1 2 1 2 ... 2 2 2 2 2 2 1 2 2
      lp               (chain, draw) float64 4.508 4.507 5.118 ... 4.63 4.667
  Attributes:
      created_at:     2024-01-23T07:45:24.002825
      arviz_version:  0.17.0

  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: (7)
    • acceptance_rate
      (chain, draw)
      float64
      0.9147 0.7716 ... 0.6239 0.9669

      array([[0.91468768, 0.7716379 , 0.84450157, ..., 0.87915957, 0.99053604,
              0.95632751],
             [0.8571711 , 0.95848208, 1.        , ..., 1.        , 0.98068892,
              0.92978242],
             [0.94753801, 0.92601269, 0.99587534, ..., 0.67847812, 0.91874785,
              0.92272868],
             [1.        , 0.99252449, 0.63301517, ..., 0.88174682, 0.62386228,
              0.96686368]])

    • step_size
      (chain, draw)
      float64
      1.183 1.183 1.183 ... 0.9726 0.9726

      array([[1.18315538, 1.18315538, 1.18315538, ..., 1.18315538, 1.18315538,
              1.18315538],
             [1.07745187, 1.07745187, 1.07745187, ..., 1.07745187, 1.07745187,
              1.07745187],
             [0.92954626, 0.92954626, 0.92954626, ..., 0.92954626, 0.92954626,
              0.92954626],
             [0.97255367, 0.97255367, 0.97255367, ..., 0.97255367, 0.97255367,
              0.97255367]])

    • 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]])

    • energy
      (chain, draw)
      float64
      5.58 5.612 5.119 ... 9.76 4.973

      array([[5.58017914, 5.61153913, 5.11860645, ..., 5.20447701, 5.1972998 ,
              4.65247233],
             [5.18502829, 5.30214299, 4.60190631, ..., 4.75489472, 4.67239329,
              5.15554611],
             [6.66427354, 5.61518822, 4.9102212 , ..., 7.16673225, 6.71951083,
              5.11791457],
             [5.93249532, 5.27518995, 9.15356324, ..., 5.44036879, 9.75997395,
              4.97279059]])

    • n_steps
      (chain, draw)
      int64
      3 3 3 1 3 1 3 1 ... 3 3 3 3 3 1 3 3

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

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

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

    • lp
      (chain, draw)
      float64
      4.508 4.507 5.118 ... 4.63 4.667

      array([[4.50795286, 4.50749236, 5.11839547, ..., 5.18281214, 4.48554271,
              4.65244764],
             [5.13045302, 4.62727449, 4.50651762, ..., 4.52430535, 4.66101964,
              4.51438521],
             [4.69756676, 4.79818193, 4.85289303, ..., 5.66177633, 4.55502287,
              4.69010656],
             [5.41221621, 4.58197845, 7.81182141, ..., 5.43805705, 4.63047395,
              4.66744693]])

  • Indexes: (2)
    • chain
      PandasIndex

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

    • draw
      PandasIndex

      PandasIndex(Index([  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: (2)

    created_at :

    2024-01-23T07:45:24.002825

    arviz_version :

    0.17.0

  
  
• observed_data

  <xarray.Dataset>
  Dimensions:    (obs_dim_0: 1)
  Coordinates:
    * obs_dim_0  (obs_dim_0) int64 0
  Data variables:
      obs        (obs_dim_0) int64 14
  Attributes:
      created_at:                 2024-01-23T07:45:24.004056
      arviz_version:              0.17.0
      inference_library:          numpyro
      inference_library_version:  0.13.2
      sampling_time:              1.29361

  xarray.Dataset

  • Dimensions:
    • obs_dim_0: 1
  • Coordinates: (1)
    • obs_dim_0
      (obs_dim_0)
      int64
      0

      array([0])

  • Data variables: (1)
    • obs
      (obs_dim_0)
      int64
      14

      array([14])

  • Indexes: (1)
    • obs_dim_0
      PandasIndex

      PandasIndex(Index([0], dtype='int64', name='obs_dim_0'))

  • Attributes: (5)

    created_at :

    2024-01-23T07:45:24.004056

    arviz_version :

    0.17.0

    inference_library :

    numpyro

    inference_library_version :

    0.13.2

    sampling_time :

    1.29361

  
  

I vari attributi dell’oggetto InferenceData possono essere estratti come si farebbe con un dict, in cui ciascuna coppia è composta da una chiave e un valore, separati dal simbolo dei due punti. In questa situazione, le chiavi sono i nomi delle variabili, mentre i valori sono rappresentati da array di Numpy. Per esempio, possiamo recuperare la traccia di campionamento dalla variabile latente theta nel modo seguente.

idata.posterior["theta"]

<xarray.DataArray 'theta' (chain: 4, draw: 1000)>
array([[0.17326516, 0.17320269, 0.20629908, ..., 0.12454636, 0.16953616,
        0.18556175],
       [0.20671803, 0.14461032, 0.17306915, ..., 0.17527697, 0.1427207 ,
        0.15351922],
       [0.14086973, 0.19343519, 0.19594644, ..., 0.1142225 , 0.17836103,
        0.14123355],
       [0.11920471, 0.14753389, 0.26728206, ..., 0.11865411, 0.14442225,
        0.14238198]])
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

xarray.DataArray

'theta'

• chain: 4
• draw: 1000

• 0.1733 0.1732 0.2063 0.1819 0.1343 ... 0.1351 0.1187 0.1444 0.1424

  array([[0.17326516, 0.17320269, 0.20629908, ..., 0.12454636, 0.16953616,
          0.18556175],
         [0.20671803, 0.14461032, 0.17306915, ..., 0.17527697, 0.1427207 ,
          0.15351922],
         [0.14086973, 0.19343519, 0.19594644, ..., 0.1142225 , 0.17836103,
          0.14123355],
         [0.11920471, 0.14753389, 0.26728206, ..., 0.11865411, 0.14442225,
          0.14238198]])

• 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])

• Indexes: (2)
  • chain
    PandasIndex

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

  • draw
    PandasIndex

    PandasIndex(Index([  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: (0)

Si noti che l’oggetto ritornato è un array di dimensioni 4 \(\times\) 1000 (sul mio computer):

idata.posterior["theta"].shape

(4, 1000)

Per visualizzare il primi 10 valori della prima catena, ad esempio, usiamo:

idata.posterior["theta"][0, 1:10]

<xarray.DataArray 'theta' (draw: 9)>
array([0.17320269, 0.20629908, 0.18193757, 0.13434879, 0.14592237,
       0.1342829 , 0.12942379, 0.2265291 , 0.23713072])
Coordinates:
    chain    int64 0
  * draw     (draw) int64 1 2 3 4 5 6 7 8 9

xarray.DataArray

'theta'

• draw: 9

• 0.1732 0.2063 0.1819 0.1343 0.1459 0.1343 0.1294 0.2265 0.2371

  array([0.17320269, 0.20629908, 0.18193757, 0.13434879, 0.14592237,
         0.1342829 , 0.12942379, 0.2265291 , 0.23713072])

• Coordinates: (2)
  • chain
    ()
    int64
    0

    array(0)

  • draw
    (draw)
    int64
    1 2 3 4 5 6 7 8 9

    array([1, 2, 3, 4, 5, 6, 7, 8, 9])

• Indexes: (1)
  • draw
    PandasIndex

    PandasIndex(Index([1, 2, 3, 4, 5, 6, 7, 8, 9], dtype='int64', name='draw'))

• Attributes: (0)

Per combinare catene e iterazioni, utilizziamo la funzione arviz.extract().

post = az.extract(idata)
post

<xarray.Dataset>
Dimensions:  (sample: 4000)
Coordinates:
  * sample   (sample) object MultiIndex
  * chain    (sample) int64 0 0 0 0 0 0 0 0 0 0 0 0 ... 3 3 3 3 3 3 3 3 3 3 3 3
  * draw     (sample) int64 0 1 2 3 4 5 6 7 ... 992 993 994 995 996 997 998 999
Data variables:
    theta    (sample) float64 0.1733 0.1732 0.2063 ... 0.1187 0.1444 0.1424
Attributes:
    created_at:     2024-01-23T07:45:24.000626
    arviz_version:  0.17.0

xarray.Dataset

• Dimensions:
  • sample: 4000
• Coordinates: (3)
  • sample
    (sample)
    object
    MultiIndex

    array([(0, 0), (0, 1), (0, 2), ..., (3, 997), (3, 998), (3, 999)], dtype=object)

  • chain
    (sample)
    int64
    0 0 0 0 0 0 0 0 ... 3 3 3 3 3 3 3 3

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

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

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

• Data variables: (1)
  • theta
    (sample)
    float64
    0.1733 0.1732 ... 0.1444 0.1424

    array([0.17326516, 0.17320269, 0.20629908, ..., 0.11865411, 0.14442225,
           0.14238198])

• Indexes: (1)
  • sample
    chain
    draw
    PandasMultiIndex

    PandasIndex(MultiIndex([(0,   0),
                (0,   1),
                (0,   2),
                (0,   3),
                (0,   4),
                (0,   5),
                (0,   6),
                (0,   7),
                (0,   8),
                (0,   9),
                ...
                (3, 990),
                (3, 991),
                (3, 992),
                (3, 993),
                (3, 994),
                (3, 995),
                (3, 996),
                (3, 997),
                (3, 998),
                (3, 999)],
               name='sample', length=4000))

• Attributes: (2)

  created_at :

  2024-01-23T07:45:24.000626

  arviz_version :

  0.17.0

La funzione extract ritorna un oggetto di classe xarray.DataArray il quale può essere trattato come un array NumPy.

post["theta"]

<xarray.DataArray 'theta' (sample: 4000)>
array([0.17326516, 0.17320269, 0.20629908, ..., 0.11865411, 0.14442225,
       0.14238198])
Coordinates:
  * sample   (sample) object MultiIndex
  * chain    (sample) int64 0 0 0 0 0 0 0 0 0 0 0 0 ... 3 3 3 3 3 3 3 3 3 3 3 3
  * draw     (sample) int64 0 1 2 3 4 5 6 7 ... 992 993 994 995 996 997 998 999

xarray.DataArray

'theta'

• sample: 4000

• 0.1733 0.1732 0.2063 0.1819 0.1343 ... 0.1351 0.1187 0.1444 0.1424

  array([0.17326516, 0.17320269, 0.20629908, ..., 0.11865411, 0.14442225,
         0.14238198])

• Coordinates: (3)
  • sample
    (sample)
    object
    MultiIndex

    array([(0, 0), (0, 1), (0, 2), ..., (3, 997), (3, 998), (3, 999)], dtype=object)

  • chain
    (sample)
    int64
    0 0 0 0 0 0 0 0 ... 3 3 3 3 3 3 3 3

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

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

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

• Indexes: (1)
  • sample
    chain
    draw
    PandasMultiIndex

    PandasIndex(MultiIndex([(0,   0),
                (0,   1),
                (0,   2),
                (0,   3),
                (0,   4),
                (0,   5),
                (0,   6),
                (0,   7),
                (0,   8),
                (0,   9),
                ...
                (3, 990),
                (3, 991),
                (3, 992),
                (3, 993),
                (3, 994),
                (3, 995),
                (3, 996),
                (3, 997),
                (3, 998),
                (3, 999)],
               name='sample', length=4000))

• Attributes: (0)

Si noti che post["theta"] contiene il numero totale (4000) dei campioni a posteriori della distribuzione di \(\theta\) escluso il burn-in:

post["theta"].shape

(4000,)

Di conseguenza, possiamo calcolare su di esso tutte le misure statistiche descrittive che si possono ottenere da un vettore di dati. Per esempio, possiamo calcolare la media a posteriori.

post.mean()

<xarray.Dataset>
Dimensions:  ()
Data variables:
    theta    float64 0.1647

xarray.Dataset

• Dimensions:
• Coordinates: (0)
• Data variables: (1)
  • theta
    ()
    float64
    0.1647

    array(0.1647053)

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

post["theta"].mean()

<xarray.DataArray 'theta' ()>
array(0.1647053)

xarray.DataArray

'theta'

• 0.1647

  array(0.1647053)

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

Possiamo calcolare la mediana a posteriori di \(\theta\).

post.median()

<xarray.Dataset>
Dimensions:  ()
Data variables:
    theta    float64 0.1626

xarray.Dataset

• Dimensions:
• Coordinates: (0)
• Data variables: (1)
  • theta
    ()
    float64
    0.1626

    array(0.16258186)

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

Oppure la deviazione standard della stima a posteriori di \(\theta\).

np.std(post["theta"])

<xarray.DataArray 'theta' ()>
array(0.03581982)

xarray.DataArray

'theta'

• 0.03582

  array(0.03581982)

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

Intervallo di credibilità

L’inferenza bayesiana tramite l’intervallo di credibilità riguarda invece la stima dell’intervallo che contiene il parametro \(\theta\) ad un dato livello di probabilità soggettiva.

Usando idata, possiamo ottenere un sommario della distribuzione a posteriori con il metodo az.summary().

az.summary(idata, hdi_prob=0.94, round_to=3)

	mean	sd	hdi_3%	hdi_97%	mcse_mean	mcse_sd	ess_bulk	ess_tail	r_hat
theta	0.165	0.036	0.099	0.233	0.001	0.001	1586.797	1643.117	1.001

Si ottiene così l’intervallo di credibilità a più alta densità a posteriori (HPD) al 94%. Questo intervallo ci informa sul fatto che, a posteriori, possiamo essere certi al 94%, che il vero valore del parametro \(\theta\) sia contenuto nell’intervallo [0.103, 0.23].

Dato che, nel caso presente, conosciamo la soluziona analitica, possiamo verificare il risultato precedente calcolando i quantili della distribuzione a posteriori Beta(18, 92) di ordine 0.03 e 0.97.

ll = stats.beta.ppf(0.03, 18, 92)
ul = stats.beta.ppf(0.97, 18, 92)
list([ll, ul])

[0.10303527075398666, 0.23457657606771784]

Verifica di ipotesi bayesiana

Un secondo tipo di inferenza bayesiana riguarda problemi in cui siamo interessati a valutare la plausibilità che il parametro \(\theta\) assuma valori contenuti in un dato intervallo di valori. Per esempio, ci potrebbe interessare l’ipotesi \(\theta > 0.5\). In questo caso, possiamo calcolare la probabilità a posteriori che \(\theta\) cada nell’intervallo di interesse, integrando la distribuzione a posteriori Beta su tale intervallo.

Nel caso dell’esempio degli artisti della Generazione X, supponiamo di essere interessati alle due seguenti ipotesi:

\[\begin{split} \begin{split} H_0: & \; \; \pi \ge 0.2 \\ H_a: & \; \; \pi < 0.2 \end{split} \end{split}\]

La nostra domanda è la seguente: Date le nostre credenze iniziali e i dati disponibili, quale importanza relativa possiamo attribuire a queste due ipotesi?

Per affrontare questa questione, iniziamo a calcolare la probabilità \(P(\theta < 0.2)\).

print(np.mean(post["theta"] < 0.2))

<xarray.DataArray 'theta' ()>
array(0.83525)

Passiamo ora a calcolare gli odds a posteriori:

post_odds = (np.mean(post["theta"] < 0.2)) / (1 - np.mean(post["theta"] < 0.2))
print(post_odds)

<xarray.DataArray 'theta' ()>
array(5.06980273)

Ciò implica che la probabilità che \(\pi\) sia inferiore al 20% è circa 6 volte superiore rispetto alla probabilità che sia al di sopra del 20%.

Questo risultato si basa solo sulle informazioni relative alla distribuzione a posteriori. Prima di avere osservato i dati del campione, avevamo una distribuzione a priori \(\operatorname{Beta}(6, 4)\), e in quel contesto avevamo una probabilità del 9% che \(H_a\) fosse vera e una probabilità del 91% che fosse falsa.

threshold = 0.2
prior_prob = stats.beta.cdf(threshold, a=alpha_prior, b=beta_prior)

prior_odds = prior_prob / (1 - prior_prob)
print(prior_odds)

0.09366320688790145

Ora possiamo combinare le informazioni degli odds a posteriori e degli odds a priori in una quantità chiamata Bayes Factor, che è semplicemente il rapporto tra le due:

BF = post_odds / prior_odds
print(BF)

<xarray.DataArray 'theta' ()>
array(54.12800714)

In conclusione, dopo aver appreso informazioni riguardo a 14 artisti appartenenti alla generazione X, le probabilità posteriori della nostra ipotesi \(H_a: \; \pi < 0.2\) sono circa 60 volte superiori rispetto alle probabilità a priori.

Questo è un esempio di test di ipotesi bayesiano.

Commenti e considerazioni finali

In questo capitolo sono state presentate le strategie per la trasformazione della distribuzione a posteriori. Sono state esplorate le diverse modalità per ottenere un intervallo di credibilità. In seguito, è stata affrontata l’analisi delle ipotesi a posteriori. Questo approccio permette di comparare in maniera relativa due ipotesi contrapposte riguardanti il parametro \(\theta\). In alcune circostanze, tale confronto viene tradotto in una misura denominata Fattore di Bayes.

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

Last updated: Tue Jan 23 2024

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

scipy     : 1.11.4
seaborn   : 0.13.0
matplotlib: 3.8.2
arviz     : 0.17.0
numpy     : 1.26.2
pandas    : 2.1.4
pymc      : 5.10.3

Watermark: 2.4.3