Posterior Predictive Checks

I controlli predittivi a posteriori (PPC) rappresentano un ottimo metodo per convalidare un modello. L’idea alla base è quella di generare dati dal modello utilizzando i parametri della distribuzione a posteriori. Questo argomento è stato trattato nel capitolo La predizione bayesiana. È importante ricordare che i PPC mirano a valutare il grado di discrepanza tra i dati generati dal modello e la vera distribuzione \(Y\). In altre parole, ciò ci permette di capire se la distribuzione a posteriori ottenuta è simile alla distribuzione teorica dei dati osservati. Nel seguente tutorial, esploreremo come utilizzare gli strumenti grafici offerti da arviz per affrontare questo problema.

Simulazione

Iniziamo a simulare dei dati.

import arviz as az
import matplotlib.pyplot as plt
import seaborn as sns
import numpy as np
import pymc as pm
from pymc import HalfNormal, Model, Normal, sample
import xarray as xr
from statsmodels.nonparametric.smoothers_lowess import lowess as sm_lowess
import warnings

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

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

Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.
Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.

Runing on PyMC v5.6.1

# 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
%matplotlib inline
%config InlineBackend.figure_format = 'retina'
%config InlineBackend.figure_format = "svg"

def standardize(series):
    """Standardize a pandas series"""
    return (series - series.mean()) / series.std()

I dati verranno simulati in modo tale da avere una relazione non lineare tra \(X\) e \(Y\).

# Size of dataset
size = 100

x = 2*np.random.random(size)
y = np.sin(x) * 12*np.exp(-x) + np.random.normal(0, 0.2, size)

zx = standardize(x)
zy = standardize(y)

sm_x, sm_y = sm_lowess(zy, zx,  frac=1./5., it=5, return_sorted=True).T
plt.plot(sm_x, sm_y, color='tomato')
plt.plot(zx, zy, 'k.')

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

È ovvio che, in un caso come questo, un modello lineare come quello che abbiamo esaminato in precedenza, non è adeguato. Ma adattiamo comunque un modello lineare ai dati proprio per verificare se i PPC saranno in grado di mettere in luce il fatto che il modello è sbagliato. Implemento qui lo stesso modello che abbiamo usato in precedenza. I dati sono standardizzati, per cui le seguenti distribuzioni a priori per i parametri sono adeguate.

with pm.Model() as model_0:
    a = pm.Normal("a", 0.0, 2.0)
    b = pm.Normal("b", 0.0, 2.0)

    mu = a + b * zx
    sigma = pm.Normal("sigma", sigma=5.0)

    pm.Normal("obs", mu=mu, sigma=sigma, observed=zy)
    idata_0 = pm.sample()

Auto-assigning NUTS sampler...

Initializing NUTS using jitter+adapt_diag...

Multiprocess sampling (4 chains in 4 jobs)

NUTS: [a, b, sigma]

Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.
Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.

Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.
Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.

Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.
Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.

Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.
Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.

100.00% [8000/8000 00:01<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 14 seconds.

Anche se il modello è sbagliato, i trace plot non mettono in evidenza alcuna anomalia di rilievo.

_ = az.plot_trace(idata_0, combined=True)
plt.tight_layout()

Generiamo ora i dati necessari per i controlli predittivi a posteriori. A tal fine, useremo una funzione PyMC dedicata a questo scopo per campionare i dati dalla distribuzione a posteriori. Questa funzione estrarrà casualmente 4000 campioni di parametri dalla traccia. Quindi, per ogni campione, estrarrà 100 numeri casuali da una distribuzione normale specificata dai valori di mu e sigma in quel campione (si veda il capitolo La predizione bayesiana.):

with model_0:
    pm.sample_posterior_predictive(
        idata_0, extend_inferencedata=True, random_seed=rng);

Sampling: [obs]

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

Ora, l’oggetto posterior_predictive in idata_0 contiene 4000 set di dati simulati, contenenti 100 campioni ciascuno, ognuno dei quali è stato calcolato usando un valore del parametro preso a caso dalla distribuzione a posteriori:

idata_0

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:
      a        (chain, draw) float64 -0.01989 -0.01834 0.03031 ... 0.1133 -0.14
      b        (chain, draw) float64 0.06934 -0.09106 0.1126 ... 0.1101 -0.09125
      sigma    (chain, draw) float64 1.058 0.9529 1.023 ... 0.9836 1.018 1.048
  Attributes:
      created_at:                 2023-09-17T09:25:16.366094
      arviz_version:              0.16.0
      inference_library:          pymc
      inference_library_version:  5.6.1
      sampling_time:              13.77347207069397
      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: (3)
    • a
      (chain, draw)
      float64
      -0.01989 -0.01834 ... 0.1133 -0.14

      array([[-0.0198888 , -0.01833798,  0.03031262, ..., -0.02325063,
              -0.01640116, -0.0438053 ],
             [ 0.04076408,  0.19006552,  0.09884467, ...,  0.02636156,
              -0.0376292 ,  0.21880037],
             [ 0.04228076, -0.02423327,  0.00754149, ...,  0.10573876,
              -0.14665813,  0.12228287],
             [-0.0133825 ,  0.01454308, -0.00669363, ...,  0.02715209,
               0.11325295, -0.14004407]])

    • b
      (chain, draw)
      float64
      0.06934 -0.09106 ... -0.09125

      array([[ 6.93391386e-02, -9.10579314e-02,  1.12614555e-01, ...,
               1.30177547e-01,  4.32077847e-03, -3.93878099e-02],
             [-1.35293304e-01,  4.44105378e-02,  7.90310564e-02, ...,
               7.60005351e-02, -1.41746949e-01,  3.67664012e-02],
             [ 1.35284992e-02,  5.68282124e-05, -2.00415124e-01, ...,
              -9.52171319e-02,  3.32626446e-02,  4.68346300e-03],
             [-3.19768922e-03, -3.90041297e-02,  2.80668952e-01, ...,
              -1.83762830e-01,  1.10120665e-01, -9.12498160e-02]])

    • sigma
      (chain, draw)
      float64
      1.058 0.9529 1.023 ... 1.018 1.048

      array([[1.05791274, 0.9529396 , 1.02267442, ..., 1.0433955 , 1.07083517,
              1.04375915],
             [1.05633956, 1.05717216, 0.98685918, ..., 0.92048997, 1.045272  ,
              1.07809577],
             [1.07866961, 0.93917479, 1.11918954, ..., 0.96336078, 1.04661687,
              0.93859473],
             [0.94249221, 1.16839729, 1.05868725, ..., 0.98357627, 1.01795652,
              1.0479325 ]])

  • 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: (6)

    created_at :

    2023-09-17T09:25:16.366094

    arviz_version :

    0.16.0

    inference_library :

    pymc

    inference_library_version :

    5.6.1

    sampling_time :

    13.77347207069397

    tuning_steps :

    1000

  
  
• posterior_predictive

  <xarray.Dataset>
  Dimensions:    (chain: 4, draw: 1000, obs_dim_2: 100)
  Coordinates:
    * chain      (chain) int64 0 1 2 3
    * draw       (draw) int64 0 1 2 3 4 5 6 7 ... 992 993 994 995 996 997 998 999
    * obs_dim_2  (obs_dim_2) int64 0 1 2 3 4 5 6 7 8 ... 92 93 94 95 96 97 98 99
  Data variables:
      obs        (chain, draw, obs_dim_2) float64 2.161 0.3026 ... 0.9169 -1.11
  Attributes:
      created_at:                 2023-09-17T09:25:17.634070
      arviz_version:              0.16.0
      inference_library:          pymc
      inference_library_version:  5.6.1

  xarray.Dataset

  • Dimensions:
    • chain: 4
    • draw: 1000
    • obs_dim_2: 100
  • 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])

    • obs_dim_2
      (obs_dim_2)
      int64
      0 1 2 3 4 5 6 ... 94 95 96 97 98 99

      array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
             18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35,
             36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53,
             54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71,
             72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89,
             90, 91, 92, 93, 94, 95, 96, 97, 98, 99])

  • Data variables: (1)
    • obs
      (chain, draw, obs_dim_2)
      float64
      2.161 0.3026 0.175 ... 0.9169 -1.11

      array([[[ 2.16126901,  0.30257532,  0.17500602, ...,  0.01884121,
               -0.87202399, -1.45221573],
              [ 0.53630361,  0.88441955,  1.5911143 , ..., -0.12150984,
               -1.46646497, -0.30585403],
              [ 1.17918428, -0.19449339, -0.12801675, ...,  0.54983531,
               -1.50736209, -0.68926855],
              ...,
              [ 1.11426912,  0.12648979, -0.2077029 , ..., -1.04616112,
                0.09995236, -1.54030372],
              [ 1.09264399,  0.73862559,  0.73473278, ..., -0.36008813,
                0.7037742 ,  1.53525618],
              [-0.77052833,  0.4546499 ,  1.00997005, ...,  0.95417078,
                0.39972887, -1.48792315]],
      
             [[ 0.67387075,  1.18344957,  0.20773195, ...,  0.88777112,
               -0.02912054, -0.46875639],
              [ 1.28332941,  0.53456538, -1.03083209, ...,  0.3566791 ,
               -1.2456953 ,  0.20533286],
              [ 1.37779415,  1.35703493,  0.30642938, ...,  0.06435585,
               -0.6634443 , -0.2644408 ],
      ...
              [ 0.83826261,  0.86829477, -0.90861022, ...,  0.5133175 ,
               -0.2418215 , -0.5326053 ],
              [ 0.73889341,  0.33492514,  0.96487244, ..., -0.50694559,
                0.57087559, -1.48475615],
              [ 0.02983323,  0.08430678,  1.15351398, ..., -0.25159083,
                2.6362519 , -0.42278967]],
      
             [[ 0.20784165, -1.96681766, -0.14493295, ..., -0.29813437,
               -0.88163185, -1.05289382],
              [-1.11682912,  2.32176197,  1.94827558, ...,  1.75279261,
                0.67419583,  0.3091305 ],
              [ 1.02596144,  1.36934456, -2.17865128, ...,  1.15409181,
                3.35421158,  0.0425367 ],
              ...,
              [-0.25703793, -0.02882394, -0.92914575, ...,  2.21514836,
               -0.73369767,  0.84158083],
              [ 1.23670177, -0.73027926, -0.56855987, ...,  0.02697896,
                1.15578681, -0.40500893],
              [ 0.04305481, -2.27353388, -1.22087098, ..., -0.83383655,
                0.91693392, -1.10963013]]])

  • Indexes: (3)
    • 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))

    • obs_dim_2
      PandasIndex

      PandasIndex(Index([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
             18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35,
             36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53,
             54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71,
             72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89,
             90, 91, 92, 93, 94, 95, 96, 97, 98, 99],
            dtype='int64', name='obs_dim_2'))

  • Attributes: (4)

    created_at :

    2023-09-17T09:25:17.634070

    arviz_version :

    0.16.0

    inference_library :

    pymc

    inference_library_version :

    5.6.1

  
  
• 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)
      reached_max_treedepth  (chain, draw) bool False False False ... False False
      index_in_trajectory    (chain, draw) int64 -2 -3 1 2 1 -2 2 ... 3 3 -1 0 2 3
      n_steps                (chain, draw) float64 3.0 3.0 3.0 3.0 ... 1.0 3.0 3.0
      largest_eigval         (chain, draw) float64 nan nan nan nan ... nan nan nan
      energy                 (chain, draw) float64 149.7 149.2 ... 151.8 150.7
      max_energy_error       (chain, draw) float64 -0.2912 0.1194 ... 0.89 0.4485
      ...                     ...
      step_size              (chain, draw) float64 0.9313 0.9313 ... 0.6795 0.6795
      lp                     (chain, draw) float64 -148.2 -148.4 ... -148.9 -149.1
      tree_depth             (chain, draw) int64 2 2 2 2 3 2 2 2 ... 2 2 2 2 1 2 2
      diverging              (chain, draw) bool False False False ... False False
      perf_counter_diff      (chain, draw) float64 0.0002626 ... 0.0002384
      smallest_eigval        (chain, draw) float64 nan nan nan nan ... nan nan nan
  Attributes:
      created_at:                 2023-09-17T09:25:16.374671
      arviz_version:              0.16.0
      inference_library:          pymc
      inference_library_version:  5.6.1
      sampling_time:              13.77347207069397
      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)
    • 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]])

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

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

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

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

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

    • energy
      (chain, draw)
      float64
      149.7 149.2 150.2 ... 151.8 150.7

      array([[149.67081431, 149.17706627, 150.22126376, ..., 149.65470377,
              149.18027942, 148.2635455 ],
             [152.95626794, 151.36152346, 149.78934772, ..., 149.45150348,
              149.39706839, 152.44148979],
             [148.49944688, 148.74717274, 151.2151111 , ..., 151.59416771,
              149.72444264, 149.55916171],
             [149.98340953, 149.97055596, 153.50079973, ..., 152.898523  ,
              151.78806311, 150.71767164]])

    • max_energy_error
      (chain, draw)
      float64
      -0.2912 0.1194 ... 0.89 0.4485

      array([[-0.29116917,  0.11944258,  0.48541121, ..., -0.07258016,
              -0.10934981, -0.096469  ],
             [ 2.35654249,  0.55605485, -0.28851645, ...,  0.21347106,
              -0.16428373,  1.24575792],
             [ 0.15391395,  0.09258276,  1.6337311 , ...,  1.68716723,
               0.12476386, -0.09755157],
             [-0.67062751,  0.56430668,  0.78002488, ...,  2.86792626,
               0.89001654,  0.44854366]])

    • perf_counter_start
      (chain, draw)
      float64
      5.216e+03 5.216e+03 ... 5.216e+03

      array([[5215.77103017, 5215.77135221, 5215.77166608, ..., 5216.09634896,
              5216.09667233, 5216.09698446],
             [5215.79256225, 5215.79348175, 5215.79384654, ..., 5216.13133662,
              5216.13164533, 5216.13195208],
             [5215.8041475 , 5215.80479683, 5215.80514396, ..., 5216.18329358,
              5216.18347508, 5216.1837865 ],
             [5216.11694442, 5216.1172695 , 5216.11759075, ..., 5216.42207   ,
              5216.42224771, 5216.42254262]])

    • energy_error
      (chain, draw)
      float64
      -0.2912 0.01941 ... 0.05964 0.2205

      array([[-0.29116917,  0.01940666,  0.11609544, ..., -0.07258016,
              -0.07278817, -0.041673  ],
             [ 0.84592438,  0.26471422, -0.28851645, ...,  0.07572827,
              -0.01198607,  0.48740737],
             [ 0.12837571, -0.07660324,  1.07746372, ...,  0.        ,
              -0.01023076, -0.03071103],
             [-0.67062751,  0.56430668,  0.26157745, ...,  0.        ,
               0.05964021,  0.22045886]])

    • process_time_diff
      (chain, draw)
      float64
      0.000263 0.000257 ... 0.000239

      array([[0.000263, 0.000257, 0.000249, ..., 0.000258, 0.000256, 0.000508],
             [0.000271, 0.000302, 0.000277, ..., 0.000253, 0.000253, 0.000253],
             [0.000267, 0.000285, 0.000268, ..., 0.000125, 0.000252, 0.000256],
             [0.000264, 0.000264, 0.000259, ..., 0.000125, 0.000241, 0.000239]])

    • step_size_bar
      (chain, draw)
      float64
      1.006 1.006 1.006 ... 1.04 1.04

      array([[1.00588644, 1.00588644, 1.00588644, ..., 1.00588644, 1.00588644,
              1.00588644],
             [1.04078549, 1.04078549, 1.04078549, ..., 1.04078549, 1.04078549,
              1.04078549],
             [1.00987425, 1.00987425, 1.00987425, ..., 1.00987425, 1.00987425,
              1.00987425],
             [1.04028776, 1.04028776, 1.04028776, ..., 1.04028776, 1.04028776,
              1.04028776]])

    • acceptance_rate
      (chain, draw)
      float64
      1.0 0.9285 0.7416 ... 0.5976 0.7632

      array([[1.        , 0.9285327 , 0.74164639, ..., 1.        , 0.96713298,
              0.99747316],
             [0.29359117, 0.78029742, 1.        , ..., 0.90073596, 1.        ,
              0.43200081],
             [0.91228955, 0.97052459, 0.43582765, ..., 0.18504297, 0.96090177,
              0.97532002],
             [0.91086228, 0.64536211, 0.61751705, ..., 0.05681663, 0.59764384,
              0.76322791]])

    • step_size
      (chain, draw)
      float64
      0.9313 0.9313 ... 0.6795 0.6795

      array([[0.93130606, 0.93130606, 0.93130606, ..., 0.93130606, 0.93130606,
              0.93130606],
             [1.12712263, 1.12712263, 1.12712263, ..., 1.12712263, 1.12712263,
              1.12712263],
             [1.14088579, 1.14088579, 1.14088579, ..., 1.14088579, 1.14088579,
              1.14088579],
             [0.67948501, 0.67948501, 0.67948501, ..., 0.67948501, 0.67948501,
              0.67948501]])

    • lp
      (chain, draw)
      float64
      -148.2 -148.4 ... -148.9 -149.1

      array([[-148.21293064, -148.37155933, -148.37858956, ..., -148.65984203,
              -148.12998704, -148.00290071],
             [-148.84453552, -149.6795031 , -148.51521239, ..., -148.78013051,
              -148.83376629, -150.33442335],
             [-148.30125517, -148.10821828, -150.44541573, ..., -148.89182097,
              -148.90760414, -148.9354167 ],
             [-148.03958083, -149.9240307 , -151.53089855, ..., -149.46386877,
              -148.91524965, -149.1488714 ]])

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

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

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

    • perf_counter_diff
      (chain, draw)
      float64
      0.0002626 0.000257 ... 0.0002384

      array([[0.00026262, 0.00025696, 0.00024862, ..., 0.00025808, 0.00025571,
              0.0005075 ],
             [0.00027017, 0.00030271, 0.00027583, ..., 0.00025333, 0.00025267,
              0.0002525 ],
             [0.00026767, 0.00028496, 0.00026804, ..., 0.00012471, 0.00025283,
              0.00025608],
             [0.00026479, 0.00026371, 0.00025838, ..., 0.00012563, 0.00024162,
              0.00023838]])

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

  • 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: (6)

    created_at :

    2023-09-17T09:25:16.374671

    arviz_version :

    0.16.0

    inference_library :

    pymc

    inference_library_version :

    5.6.1

    sampling_time :

    13.77347207069397

    tuning_steps :

    1000

  
  
• observed_data

  <xarray.Dataset>
  Dimensions:    (obs_dim_0: 100)
  Coordinates:
    * obs_dim_0  (obs_dim_0) int64 0 1 2 3 4 5 6 7 8 ... 92 93 94 95 96 97 98 99
  Data variables:
      obs        (obs_dim_0) float64 -0.3292 0.4671 -0.9968 ... -0.5284 0.7708
  Attributes:
      created_at:                 2023-09-17T09:25:16.377642
      arviz_version:              0.16.0
      inference_library:          pymc
      inference_library_version:  5.6.1

  xarray.Dataset

  • Dimensions:
    • obs_dim_0: 100
  • Coordinates: (1)
    • obs_dim_0
      (obs_dim_0)
      int64
      0 1 2 3 4 5 6 ... 94 95 96 97 98 99

      array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
             18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35,
             36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53,
             54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71,
             72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89,
             90, 91, 92, 93, 94, 95, 96, 97, 98, 99])

  • Data variables: (1)
    • obs
      (obs_dim_0)
      float64
      -0.3292 0.4671 ... -0.5284 0.7708

      array([-0.3291907 ,  0.46708713, -0.996844  ,  0.83664244,  0.67989973,
              0.78859201, -0.60524466,  1.30932419,  0.64194125, -1.83574714,
             -0.78731874, -2.24216547, -1.76335824,  1.13966211,  1.31558497,
              1.15990829, -1.14220628,  0.0785071 ,  1.30449163, -0.22133161,
             -0.68852686, -1.16099624, -0.37830463,  1.25385287, -0.60879757,
             -1.12491657, -0.1169575 ,  1.34167288,  0.74128552, -1.14185325,
              0.32598647,  0.68042776,  0.43983195, -0.71900203, -0.88389458,
             -0.99306937, -1.58880483, -1.19912518,  1.15950733, -0.4343844 ,
             -0.6622863 ,  0.59881534, -0.91460161,  0.20228693, -2.31111612,
              1.09005131, -2.37265136, -1.38133487,  1.0347458 ,  0.65866406,
              0.34366857,  0.81854309,  1.41207653, -1.97814241, -0.55824128,
              1.04900245,  1.24897761,  0.87272623, -0.41135846,  0.3973573 ,
             -0.03014266, -0.51618908,  0.2504479 ,  0.23626335,  0.75723634,
             -1.22082725,  0.31655698,  1.10262968,  0.9933467 ,  1.026237  ,
              0.89216137, -0.54836509, -0.16105795, -1.31113892,  0.91226407,
              0.07336096,  1.04258891, -0.15748385, -0.5453532 ,  0.91895186,
              1.1386864 ,  1.13750579, -0.65059937, -0.31561081, -1.30861348,
              1.17378748,  1.09887035, -0.61808358, -0.49020392,  0.93111712,
              0.82091202, -0.30062435, -1.68255229,  1.02331241, -0.02146963,
              1.28584927, -1.2819217 , -0.0535586 , -0.52843367,  0.77079488])

  • Indexes: (1)
    • obs_dim_0
      PandasIndex

      PandasIndex(Index([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
             18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35,
             36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53,
             54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71,
             72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89,
             90, 91, 92, 93, 94, 95, 96, 97, 98, 99],
            dtype='int64', name='obs_dim_0'))

  • Attributes: (4)

    created_at :

    2023-09-17T09:25:16.377642

    arviz_version :

    0.16.0

    inference_library :

    pymc

    inference_library_version :

    5.6.1

  
  

A questo punto possiamo utilizzare la funzione az.plot_ppc() per determinare se il modello è almeno in grado di riprodurre i dati osservati:

_ = az.plot_ppc(idata_0, num_pp_samples=100)

Il PPC plot mostra come vi sia pochissima corrispondenza tra la distribuzione dei dati osservati (in nero) e la distribuzione dei dati predetti dal modello. Questo era atteso, dato che abbiamo adattato un modello del tutto sbagliato per i dati a disposizione.

Proseguiamo con la simulazione. I dati sono stati generati con una funzione non lineare. In generale, quella funzione è sconosciuta. Qualunque funzione non lineare, però, può essere approssimata da un modello polinomiale. Un modello di regressione polinomiale si ottiene inserendo nel modello dei nuovi predittori, ciascuno ottenuto elevando la variabile originale a una potenza:

\[ \mathbb{E}(Y) = \alpha + \beta_1 x + \beta_2 x^2 + \beta_3 x^3 + \beta_4 x^4 + \dots \]

Questo modello, anche se non descrive una relazione lineare tra \(X\) e \(Y\), è funzione lineare dei coefficienti ignoti \(\alpha, \beta_1, \beta_2, \dots\). Il grado del polinomio dipende dal tipo di curva che vogliamo approssimare. Proviamo qui con un modello polinomiale di ordine 4.

zx2 = zx**2
zx3 = zx**3
zx4 = zx**4

with pm.Model() as model_1:
    a = pm.Normal("a", 0.0, 3.0)
    b1 = pm.Normal("b1", 0.0, 3.0)
    b2 = pm.Normal("b2", 0.0, 3.0)
    b3 = pm.Normal("b3", 0.0, 3.0)
    b4 = pm.Normal("b4", 0.0, 3.0)
    
    mu = a + b1 * zx + b2 * zx2 + b3 * zx3 + b4 * zx4
    sigma = pm.HalfNormal("sigma", sigma=5.0)

    pm.Normal("obs", mu=mu, sigma=sigma, observed=zy)
    idata_1 = pm.sample()

Auto-assigning NUTS sampler...

Initializing NUTS using jitter+adapt_diag...

Multiprocess sampling (4 chains in 4 jobs)

NUTS: [a, b1, b2, b3, b4, sigma]

Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.
Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.

Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.
Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.

Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.
Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.

Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.
Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.

100.00% [8000/8000 00:04<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 18 seconds.

_ = az.plot_trace(idata_1, combined=True)
plt.tight_layout()

az.summary(idata_1)

	mean	sd	hdi_3%	hdi_97%	mcse_mean	mcse_sd	ess_bulk	ess_tail	r_hat
a	1.081	0.038	1.005	1.146	0.001	0.001	2190.0	2399.0	1.0
b1	-0.748	0.056	-0.850	-0.643	0.001	0.001	1663.0	1739.0	1.0
b2	-0.912	0.081	-1.062	-0.763	0.002	0.001	1558.0	2009.0	1.0
b3	0.473	0.031	0.416	0.534	0.001	0.001	1727.0	2059.0	1.0
b4	-0.111	0.034	-0.176	-0.049	0.001	0.001	1632.0	2205.0	1.0
sigma	0.202	0.015	0.176	0.231	0.000	0.000	3009.0	2558.0	1.0

Anche in questo caso, campioniamo dalla distribuzione a posteriori dei parametri e simuliamo i dati possibili futuri. L’esame del PPC plot, in questo secondo caso, mostra una buona corrispondenza tra le predizioni del modello e i dati osservati.

with model_1:
    pm.sample_posterior_predictive(
        idata_1, extend_inferencedata=True, random_seed=rng)

_ = az.plot_ppc(idata_1, num_pp_samples=100)

Sampling: [obs]

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

Si noti che la curva nera rappresenta la distribuzione della variabile dipendente, come risulta anche dal seguente KDE plot.

sns.kdeplot(zy)

<Axes: ylabel='Density'>

Poniamoci ora il problema di generare un grafico che mostra la relazione prevista tra il predittore (\(x\)) e i valori della \(y\) previsti dal modello.

with model_1:
    pm.sample_posterior_predictive(
        idata_1, extend_inferencedata=True, random_seed=rng)

post = idata_1.posterior
mu_pp = post["a"] + post["b1"] * xr.DataArray(zx, dims=["obs_id"]) + \
    post["b2"] * xr.DataArray(zx2, dims=["obs_id"]) + \
    post["b3"] * xr.DataArray(zx3, dims=["obs_id"]) + \
    post["b4"] * xr.DataArray(zx4, dims=["obs_id"]) 

_, ax = plt.subplots()

ax.plot(
    zx, mu_pp.mean(("chain", "draw")), '.', label="Mean outcome", color="C1", alpha=0.6
)
ax.scatter(zx, idata_1.observed_data["obs"])
az.plot_hdi(zx, idata_1.posterior_predictive["obs"])

ax.set_xlabel("Predictor (stdz)")
ax.set_ylabel("Outcome (stdz)")

Sampling: [obs]

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

Text(0, 0.5, 'Outcome (stdz)')

per questo secondo modello vediamo come i dati predetti (la banda evidenziata) approssimano da vicino la relazione tra i dati osservati e la \(x\).

Ripetiamo questo procedimento usando i dati del modello sbagliato.

with model_0:
    pm.sample_posterior_predictive(
        idata_0, extend_inferencedata=True, random_seed=rng)

post = idata_0.posterior
mu_pp = post["a"] + post["b"] * xr.DataArray(zx, dims=["obs_id"]) 

_, ax = plt.subplots()

ax.plot(
    zx, mu_pp.mean(("chain", "draw")), '.', label="Mean outcome", color="C1", alpha=0.6
)
ax.scatter(zx, idata_0.observed_data["obs"])
az.plot_hdi(zx, idata_0.posterior_predictive["obs"])

ax.set_xlabel("Predictor (stdz)")
ax.set_ylabel("Outcome (stdz)")

Sampling: [obs]

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

Text(0, 0.5, 'Outcome (stdz)')

Si noti come, in questo caso, vi è un’enorme incertezza nella predizone del modello. E, ovviamente, i dati predetti non rendono conto, in nessun modo, della relazione tra i valori osservati \(x, y\).

In conlusione, dunque, i PPC plot sono un’utile strumento per il confronto tra modelli.

Watermark

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

Last updated: Sun Sep 17 2023

Python implementation: CPython
Python version       : 3.11.5
IPython version      : 8.15.0

numpy     : 1.25.2
pymc      : 5.6.1
matplotlib: 3.7.2
arviz     : 0.16.0
xarray    : 2023.6.0
seaborn   : 0.12.2

Watermark: 2.4.3