Use a standard dataset for this example

In [1]:

import arviz as az
import numpy as np  # For vectorized math operations
import pandas as pd  # For file input/output
import pymc as pm
import pytensor.tensor as pt

from matplotlib import pyplot as plt
from matplotlib.lines import Line2D

WARNING (pytensor.tensor.blas): Using NumPy C-API based implementation for BLAS functions.

In [2]:

%config InlineBackend.figure_format = 'retina'  # high resolution figures
az.style.use("arviz-darkgrid")
plt.rcParams["font.family"] = "Latin Modern Roman"

rng = np.random.default_rng(42)

In [3]:

try:
    wide_heating_df = pd.read_csv("../data/heating_data_r.csv")
except:
    wide_heating_df = pd.read_csv(pm.get_data("heating_data_r.csv"))

wide_heating_df[wide_heating_df["idcase"] == 1]

	idcase	depvar	ic.gc	ic.gr	ic.ec	ic.er	ic.hp	oc.gc	oc.gr	oc.ec	oc.er	oc.hp	income	agehed	rooms	region
0	1	gc	866.0	962.64	859.9	995.76	1135.5	199.69	151.72	553.34	505.6	237.88	7	25	6	ncostl

In [4]:

# pivoted data that is more amenable to analysis
 
try:
    long_heating_df = pd.read_csv("../data/long_heating_data.csv")
except:
    long_heating_df = pd.read_csv(pm.get_data("long_heating_data.csv"))

columns = [c for c in long_heating_df.columns if c != "Unnamed: 0"]
long_heating_df[long_heating_df["idcase"] == 1][columns]

	idcase	alt_id	choice	depvar	income	agehed	rooms	region	installation_costs	operating_costs
0	1	1	1	gc	7	25	6	ncostl	866.00	199.69
1	1	2	0	gc	7	25	6	ncostl	962.64	151.72
2	1	3	0	gc	7	25	6	ncostl	859.90	553.34
3	1	4	0	gc	7	25	6	ncostl	995.76	505.60
4	1	5	0	gc	7	25	6	ncostl	1135.50	237.88

In [5]:

N = wide_heating_df.shape[0]
observed = pd.Categorical(wide_heating_df["depvar"]).codes
coords = {
    "alts_probs": ["ec", "er", "gc", "gr", "hp"],
    "obs": range(N),
}

with pm.Model(coords=coords) as model_1:
    beta_ic = pm.Normal("beta_ic", 0, 1)
    beta_oc = pm.Normal("beta_oc", 0, 1)

    ## Construct Utility matrix and Pivot
    u0 = beta_ic * wide_heating_df["ic.ec"] + beta_oc * wide_heating_df["oc.ec"]
    u1 = beta_ic * wide_heating_df["ic.er"] + beta_oc * wide_heating_df["oc.er"]
    u2 = beta_ic * wide_heating_df["ic.gc"] + beta_oc * wide_heating_df["oc.gc"]
    u3 = beta_ic * wide_heating_df["ic.gr"] + beta_oc * wide_heating_df["oc.gr"]
    u4 = np.zeros(N)  # Outside Good
    s = pm.math.stack([u0, u1, u2, u3, u4]).T

    ## Apply Softmax Transform
    p_ = pm.Deterministic("p", pm.math.softmax(s, axis=1), dims=("obs", "alts_probs"))

    ## Likelihood
    choice_obs = pm.Categorical("y_cat", p=p_, observed=observed, dims="obs")

    idata_m1 = pm.sample_prior_predictive()
    idata_m1.extend(
        pm.sample(nuts_sampler="numpyro", idata_kwargs={"log_likelihood": True}, random_seed=101)
    )
    idata_m1.extend(pm.sample_posterior_predictive(idata_m1))

pm.model_to_graphviz(model_1)

Sampling: [beta_ic, beta_oc, y_cat]

Sampling: [y_cat]

In [6]:

idata_m1

arviz.InferenceData

• posterior

  <xarray.Dataset> Size: 144MB
  Dimensions:     (chain: 4, draw: 1000, obs: 900, alts_probs: 5)
  Coordinates:
    * chain       (chain) int64 32B 0 1 2 3
    * draw        (draw) int64 8kB 0 1 2 3 4 5 6 7 ... 993 994 995 996 997 998 999
    * obs         (obs) int64 7kB 0 1 2 3 4 5 6 7 ... 893 894 895 896 897 898 899
    * alts_probs  (alts_probs) <U2 40B 'ec' 'er' 'gc' 'gr' 'hp'
  Data variables:
      beta_ic     (chain, draw) float64 32kB 0.001832 0.001817 ... 0.001948
      beta_oc     (chain, draw) float64 32kB -0.003705 -0.003709 ... -0.004581
      p           (chain, draw, obs, alts_probs) float64 144MB 0.07556 ... 0.1102
  Attributes:
      created_at:     2024-06-23T20:32:15.207517+00:00
      arviz_version:  0.18.0

  xarray.Dataset

  • Dimensions:
    • chain: 4
    • draw: 1000
    • obs: 900
    • alts_probs: 5
  • Coordinates: (4)
    • 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
      (obs)
      int64
      0 1 2 3 4 5 ... 895 896 897 898 899

      array([  0,   1,   2, ..., 897, 898, 899])

    • alts_probs
      (alts_probs)
      <U2
      'ec' 'er' 'gc' 'gr' 'hp'

      array(['ec', 'er', 'gc', 'gr', 'hp'], dtype='<U2')

  • Data variables: (3)
    • beta_ic
      (chain, draw)
      float64
      0.001832 0.001817 ... 0.001948

      array([[0.00183239, 0.00181733, 0.00180854, ..., 0.00181192, 0.00175778,
              0.00210304],
             [0.00210654, 0.00210654, 0.00197494, ..., 0.0019049 , 0.00187513,
              0.00181919],
             [0.00181312, 0.00179163, 0.00173467, ..., 0.002031  , 0.00212423,
              0.00171305],
             [0.00185484, 0.00177189, 0.00187684, ..., 0.00206325, 0.00207469,
              0.00194822]])

    • beta_oc
      (chain, draw)
      float64
      -0.003705 -0.003709 ... -0.004581

      array([[-0.00370539, -0.00370935, -0.00431479, ..., -0.00461288,
              -0.00460595, -0.00406624],
             [-0.00443282, -0.00443282, -0.00414015, ..., -0.00426425,
              -0.00412309, -0.0040975 ],
             [-0.00395141, -0.00370944, -0.00377519, ..., -0.0049641 ,
              -0.00476948, -0.0036491 ],
             [-0.00381868, -0.00404807, -0.00393956, ..., -0.0042988 ,
              -0.00462425, -0.00458066]])

    • p
      (chain, draw, obs, alts_probs)
      float64
      0.07556 0.1157 ... 0.428 0.1102

      array([[[[0.07556231, 0.115677  , 0.28331138, 0.4039845 , 0.1214648 ],
               [0.09409705, 0.12757696, 0.30515943, 0.32297186, 0.15019469],
               [0.10440213, 0.16496078, 0.23067664, 0.3579255 , 0.14203495],
               ...,
               [0.08474306, 0.13638204, 0.23587751, 0.43633533, 0.10666206],
               [0.07042229, 0.18021517, 0.22561794, 0.42184934, 0.10189526],
               [0.08717172, 0.13668048, 0.27121101, 0.40382201, 0.10111478]],
      
              [[0.07540941, 0.11522867, 0.28310845, 0.40318444, 0.12306902],
               [0.09378437, 0.12698266, 0.30488553, 0.32253151, 0.15181593],
               [0.1041803 , 0.16418621, 0.23085429, 0.35725122, 0.14352798],
               ...,
               [0.08467137, 0.13587604, 0.23615981, 0.43515429, 0.10813849],
               [0.07046016, 0.17950217, 0.22578659, 0.42086267, 0.10338842],
               [0.08701837, 0.13621822, 0.27140404, 0.40272926, 0.10263011]],
      
              [[0.06122794, 0.0961875 , 0.28473604, 0.41709752, 0.140751  ],
               [0.07754801, 0.10707428, 0.3120934 , 0.3300667 , 0.17321761],
               [0.08971942, 0.14413631, 0.23485909, 0.36901728, 0.1622679 ],
               ...,
      ...
               ...,
               [0.07148748, 0.12262564, 0.23524756, 0.47222741, 0.09841191],
               [0.05809065, 0.1693353 , 0.22340738, 0.45541326, 0.09375342],
               [0.07380963, 0.12341259, 0.27669202, 0.43293536, 0.0931504 ]],
      
              [[0.0557185 , 0.09210605, 0.28955146, 0.44171102, 0.12091297],
               [0.07321771, 0.10485256, 0.32256519, 0.34396418, 0.15540036],
               [0.08500087, 0.14479629, 0.23452216, 0.39027898, 0.14540171],
               ...,
               [0.06514664, 0.11311106, 0.2365109 , 0.48142957, 0.10380182],
               [0.05287363, 0.15948594, 0.22320083, 0.4650963 , 0.09934329],
               [0.0669652 , 0.11434375, 0.28043246, 0.43930699, 0.09895161]],
      
              [[0.05622438, 0.09116917, 0.28748885, 0.43233235, 0.13278524],
               [0.07282079, 0.10284997, 0.31870176, 0.33851648, 0.16711101],
               [0.08500239, 0.14132398, 0.2353021 , 0.38212832, 0.15624321],
               ...,
               [0.06612004, 0.11180751, 0.23856091, 0.46873109, 0.11478046],
               [0.05429461, 0.15656816, 0.22464392, 0.4540332 , 0.11046011],
               [0.06750051, 0.11316997, 0.28109683, 0.4280232 , 0.11020948]]]])

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

      PandasIndex(Index([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,
             ...
             890, 891, 892, 893, 894, 895, 896, 897, 898, 899],
            dtype='int64', name='obs', length=900))

    • alts_probs
      PandasIndex

      PandasIndex(Index(['ec', 'er', 'gc', 'gr', 'hp'], dtype='object', name='alts_probs'))

  • Attributes: (2)

    created_at :

    2024-06-23T20:32:15.207517+00:00

    arviz_version :

    0.18.0

  
  
• posterior_predictive

  <xarray.Dataset> Size: 29MB
  Dimensions:  (chain: 4, draw: 1000, obs: 900)
  Coordinates:
    * chain    (chain) int64 32B 0 1 2 3
    * draw     (draw) int64 8kB 0 1 2 3 4 5 6 7 ... 993 994 995 996 997 998 999
    * obs      (obs) int64 7kB 0 1 2 3 4 5 6 7 ... 892 893 894 895 896 897 898 899
  Data variables:
      y_cat    (chain, draw, obs) int64 29MB 3 3 2 3 3 4 3 2 2 ... 3 3 3 2 3 3 4 1
  Attributes:
      created_at:                 2024-06-23T20:32:21.823805+00:00
      arviz_version:              0.18.0
      inference_library:          pymc
      inference_library_version:  5.15.1

  xarray.Dataset

  • Dimensions:
    • chain: 4
    • draw: 1000
    • obs: 900
  • 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
      (obs)
      int64
      0 1 2 3 4 5 ... 895 896 897 898 899

      array([  0,   1,   2, ..., 897, 898, 899])

  • Data variables: (1)
    • y_cat
      (chain, draw, obs)
      int64
      3 3 2 3 3 4 3 2 ... 3 3 3 2 3 3 4 1

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

  • 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
      PandasIndex

      PandasIndex(Index([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,
             ...
             890, 891, 892, 893, 894, 895, 896, 897, 898, 899],
            dtype='int64', name='obs', length=900))

  • Attributes: (4)

    created_at :

    2024-06-23T20:32:21.823805+00:00

    arviz_version :

    0.18.0

    inference_library :

    pymc

    inference_library_version :

    5.15.1

  
  
• log_likelihood

  <xarray.Dataset> Size: 29MB
  Dimensions:      (chain: 4, draw: 1000, y_cat_dim_0: 1, y_cat_dim_1: 900)
  Coordinates:
    * chain        (chain) int64 32B 0 1 2 3
    * draw         (draw) int64 8kB 0 1 2 3 4 5 6 ... 993 994 995 996 997 998 999
    * y_cat_dim_0  (y_cat_dim_0) int64 8B 0
    * y_cat_dim_1  (y_cat_dim_1) int64 7kB 0 1 2 3 4 5 ... 894 895 896 897 898 899
  Data variables:
      y_cat        (chain, draw, y_cat_dim_0, y_cat_dim_1) float64 29MB -1.261 ...
  Attributes:
      created_at:     2024-06-23T20:32:15.210389+00:00
      arviz_version:  0.18.0

  xarray.Dataset

  • Dimensions:
    • chain: 4
    • draw: 1000
    • y_cat_dim_0: 1
    • y_cat_dim_1: 900
  • Coordinates: (4)
    • 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])

    • y_cat_dim_0
      (y_cat_dim_0)
      int64
      0

      array([0])

    • y_cat_dim_1
      (y_cat_dim_1)
      int64
      0 1 2 3 4 5 ... 895 896 897 898 899

      array([  0,   1,   2, ..., 897, 898, 899])

  • Data variables: (1)
    • y_cat
      (chain, draw, y_cat_dim_0, y_cat_dim_1)
      float64
      -1.261 -1.187 ... -1.493 -1.269

      array([[[[-1.2612087 , -1.18692092, -1.46673836, ..., -1.44444265,
                -1.48891222, -1.30485813]],
      
              [[-1.26192523, -1.18781889, -1.46596855, ..., -1.44324653,
                -1.48816503, -1.30414666]],
      
              [[-1.25619272, -1.16445278, -1.44876958, ..., -1.42865775,
                -1.48677005, -1.27635524]],
      
              ...,
      
              [[-1.25606194, -1.15560848, -1.44314675, ..., -1.42458852,
                -1.48873712, -1.26557903]],
      
              [[-1.26036161, -1.16141909, -1.44246931, ..., -1.4220514 ,
                -1.48744086, -1.26553536]],
      
              [[-1.2432454 , -1.15168259, -1.46743232, ..., -1.45639165,
                -1.50228702, -1.29804564]]],
      
      ...
      
             [[[-1.25808289, -1.1799529 , -1.46386658, ..., -1.4428013 ,
                -1.4892258 , -1.2997839 ]],
      
              [[-1.26018275, -1.17702843, -1.45418793, ..., -1.43166875,
                -1.48506219, -1.2867652 ]],
      
              [[-1.25512316, -1.17290448, -1.4609046 , ..., -1.44114666,
                -1.48966503, -1.29452733]],
      
              ...,
      
              [[-1.24190813, -1.14419259, -1.45813717, ..., -1.44711688,
                -1.49875835, -1.28485022]],
      
              [[-1.23942224, -1.13145002, -1.45020518, ..., -1.441761  ,
                -1.49968331, -1.27142236]],
      
              [[-1.2465712 , -1.14349955, -1.44688504, ..., -1.43313061,
                -1.49323869, -1.26905606]]]])

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

    • y_cat_dim_0
      PandasIndex

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

    • y_cat_dim_1
      PandasIndex

      PandasIndex(Index([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,
             ...
             890, 891, 892, 893, 894, 895, 896, 897, 898, 899],
            dtype='int64', name='y_cat_dim_1', length=900))

  • Attributes: (2)

    created_at :

    2024-06-23T20:32:15.210389+00:00

    arviz_version :

    0.18.0

  
  
• sample_stats

  <xarray.Dataset> Size: 204kB
  Dimensions:          (chain: 4, draw: 1000)
  Coordinates:
    * chain            (chain) int64 32B 0 1 2 3
    * draw             (draw) int64 8kB 0 1 2 3 4 5 6 ... 994 995 996 997 998 999
  Data variables:
      acceptance_rate  (chain, draw) float64 32kB 0.9628 1.0 ... 0.9096 0.8808
      diverging        (chain, draw) bool 4kB False False False ... False False
      energy           (chain, draw) float64 32kB 1.312e+03 ... 1.312e+03
      lp               (chain, draw) float64 32kB 1.311e+03 ... 1.311e+03
      n_steps          (chain, draw) int64 32kB 3 1 7 1 3 7 3 7 ... 7 3 3 3 3 7 3
      step_size        (chain, draw) float64 32kB 0.03082 0.03082 ... 0.0351
      tree_depth       (chain, draw) int64 32kB 2 1 3 1 2 3 2 3 ... 3 2 2 2 2 3 2
  Attributes:
      created_at:     2024-06-23T20:32:15.209472+00:00
      arviz_version:  0.18.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.9628 1.0 0.9231 ... 0.9096 0.8808

      array([[0.96283342, 1.        , 0.9230903 , ..., 0.75914853, 0.76587015,
              0.9716857 ],
             [0.98301043, 0.64214544, 0.71499472, ..., 0.99954378, 0.98732408,
              0.93728333],
             [0.97384869, 0.5090271 , 0.98637995, ..., 0.97621634, 1.        ,
              1.        ],
             [0.95595037, 0.9920607 , 0.95920415, ..., 0.88600583, 0.90960034,
              0.88083658]])

    • 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
      1.312e+03 1.311e+03 ... 1.312e+03

      array([[1311.62704271, 1311.40908264, 1312.02358383, ..., 1312.94174785,
              1313.8621048 , 1314.56025091],
             [1311.40789762, 1312.88562855, 1313.431072  , ..., 1309.88449978,
              1310.48997299, 1310.14581203],
             [1310.74492576, 1315.67927381, 1311.61594878, ..., 1313.55932677,
              1313.28207282, 1311.9541163 ],
             [1310.76018534, 1310.87052439, 1310.82664292, ..., 1312.72419099,
              1312.61614047, 1312.07832766]])

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

      array([[1311.37579102, 1311.27521994, 1310.55444545, ..., 1312.52991085,
              1313.54709578, 1312.22538983],
             [1311.04717019, 1311.04717019, 1310.18177297, ..., 1309.86136009,
              1309.84522774, 1309.9889148 ],
             [1310.16298245, 1311.19549703, 1310.93490151, ..., 1313.3683328 ,
              1311.97755372, 1311.60368489],
             [1310.75790973, 1310.26335757, 1310.27104901, ..., 1310.67819214,
              1311.08020188, 1310.81157147]])

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

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

    • step_size
      (chain, draw)
      float64
      0.03082 0.03082 ... 0.0351 0.0351

      array([[0.03081926, 0.03081926, 0.03081926, ..., 0.03081926, 0.03081926,
              0.03081926],
             [0.0387321 , 0.0387321 , 0.0387321 , ..., 0.0387321 , 0.0387321 ,
              0.0387321 ],
             [0.03365488, 0.03365488, 0.03365488, ..., 0.03365488, 0.03365488,
              0.03365488],
             [0.03509999, 0.03509999, 0.03509999, ..., 0.03509999, 0.03509999,
              0.03509999]])

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

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

  • 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-06-23T20:32:15.209472+00:00

    arviz_version :

    0.18.0

  
  
• prior

  <xarray.Dataset> Size: 18MB
  Dimensions:     (chain: 1, draw: 500, obs: 900, alts_probs: 5)
  Coordinates:
    * chain       (chain) int64 8B 0
    * draw        (draw) int64 4kB 0 1 2 3 4 5 6 7 ... 493 494 495 496 497 498 499
    * obs         (obs) int64 7kB 0 1 2 3 4 5 6 7 ... 893 894 895 896 897 898 899
    * alts_probs  (alts_probs) <U2 40B 'ec' 'er' 'gc' 'gr' 'hp'
  Data variables:
      beta_ic     (chain, draw) float64 4kB 1.383 -0.07456 ... -0.921 1.762
      beta_oc     (chain, draw) float64 4kB -0.2187 1.131 ... -1.086 -0.03639
      p           (chain, draw, obs, alts_probs) float64 18MB 1.347e-100 ... 0.0
  Attributes:
      created_at:                 2024-06-23T20:32:10.756506+00:00
      arviz_version:              0.18.0
      inference_library:          pymc
      inference_library_version:  5.15.1

  xarray.Dataset

  • Dimensions:
    • chain: 1
    • draw: 500
    • obs: 900
    • alts_probs: 5
  • Coordinates: (4)
    • chain
      (chain)
      int64
      0

      array([0])

    • draw
      (draw)
      int64
      0 1 2 3 4 5 ... 495 496 497 498 499

      array([  0,   1,   2, ..., 497, 498, 499])

    • obs
      (obs)
      int64
      0 1 2 3 4 5 ... 895 896 897 898 899

      array([  0,   1,   2, ..., 897, 898, 899])

    • alts_probs
      (alts_probs)
      <U2
      'ec' 'er' 'gc' 'gr' 'hp'

      array(['ec', 'er', 'gc', 'gr', 'hp'], dtype='<U2')

  • Data variables: (3)
    • beta_ic
      (chain, draw)
      float64
      1.383 -0.07456 ... -0.921 1.762

      array([[ 1.38326650e+00, -7.45587282e-02, -5.81792481e-01,
               1.64583486e+00,  9.81388918e-01,  3.20754793e-01,
              -1.44680278e+00, -2.06025789e-01,  8.64762917e-01,
               1.73868916e+00,  3.60089734e-01,  3.88279399e-01,
               4.46724859e-01, -1.92723685e-01,  3.77030938e-02,
               2.26743914e+00, -1.66789474e-01,  1.26919707e+00,
              -5.12032636e-01, -5.55252577e-01, -4.96782131e-02,
              -2.52394339e-01,  1.42892434e+00,  3.51714532e-01,
              -2.22922874e-01, -3.50169472e-01,  1.97139329e+00,
               7.95475057e-01,  6.25534373e-01,  4.41159488e-01,
              -8.67342839e-01,  5.07451471e-01, -1.94339389e-01,
               7.84437654e-01,  9.00848272e-01, -1.58365931e+00,
              -3.94153497e-01,  4.73934137e-01, -2.32096110e-01,
               1.52679971e+00, -7.33090872e-01,  5.40440435e-02,
               4.20206750e-01,  1.68133009e-01,  8.98569623e-01,
              -1.20163327e+00, -9.75584243e-01, -1.18505918e+00,
               3.22548915e-01, -6.53172821e-01, -2.50655596e+00,
               6.38678217e-01, -1.14182974e+00,  1.16539989e+00,
               1.63830149e+00,  6.32884908e-01, -2.42133428e-01,
               2.01297425e-01, -4.94791176e-02,  3.85060974e-02,
      ...
              -1.56106162e+00, -2.84616072e-01,  1.08531916e+00,
               4.81523002e-01,  7.11481104e-01, -5.44618716e-01,
              -1.25751968e+00, -6.89922007e-01, -1.68533020e-01,
               1.14727326e-01,  8.37509204e-01,  3.64710178e-01,
               5.81156618e-02, -5.94463594e-02, -4.25593140e-03,
              -5.75202823e-01, -1.94502686e+00, -7.01523024e-01,
               3.64551865e-01, -4.26902836e-01, -8.22513091e-02,
               1.26542226e+00, -2.29849186e+00,  2.77434019e-01,
               5.18548735e-01, -8.32959448e-01,  9.90450125e-01,
               1.65323866e+00, -4.51926171e-01,  1.24054403e+00,
               1.01275885e+00, -1.42925187e+00, -1.12880836e-01,
              -3.48517306e-01, -8.73964461e-01, -4.05496974e-01,
               1.25652380e+00, -2.11422792e-02,  9.51226150e-01,
               1.19009122e+00, -1.09275551e-01,  8.38360828e-01,
               1.77667865e-01, -1.11184601e+00, -2.60523819e-01,
               9.41971128e-01,  1.34351006e+00,  2.35053486e-03,
              -1.24979756e+00, -1.46555611e+00, -9.42550642e-01,
               3.34928240e-01, -1.76595489e+00,  1.86342911e-01,
               2.20608550e-01, -2.34437706e-01,  1.11062114e+00,
              -9.21044778e-01,  1.76173066e+00]])

    • beta_oc
      (chain, draw)
      float64
      -0.2187 1.131 ... -1.086 -0.03639

      array([[-2.18724272e-01,  1.13098028e+00, -9.37282697e-01,
               2.64177807e-01, -1.17734793e+00, -1.37330362e+00,
              -1.73970760e-01, -1.36896541e+00,  3.85066116e-02,
               8.64166726e-01,  1.68653751e+00, -5.45321315e-01,
               1.01895759e+00,  8.90660529e-02,  4.44408281e-01,
               1.96294788e-01,  6.89205806e-01,  1.14307422e+00,
               3.35997470e-01, -1.61744354e+00,  9.89014201e-01,
              -2.36364200e+00, -7.47573088e-01,  1.07466100e+00,
              -2.17134798e+00,  1.29394173e-01, -5.31345888e-01,
               1.02857117e-01,  1.41281542e+00,  2.17402322e+00,
               7.17895189e-01, -7.74930820e-01,  1.24237597e+00,
              -4.65215265e-01,  1.37410741e+00,  8.68770005e-01,
              -2.38616653e-01,  5.33523384e-01, -8.05903979e-01,
              -6.98977110e-01, -1.44637213e+00, -1.09441735e-01,
              -5.03454300e-01,  6.20235572e-01,  4.13557688e-01,
               1.06777670e+00, -5.84416370e-01,  2.17311742e+00,
              -1.42837875e+00,  5.28116978e-01, -4.96708329e-01,
              -6.29428281e-01, -6.06128256e-01,  7.84153504e-01,
               1.53397302e-01, -1.02352567e+00, -3.43940374e-01,
              -3.70735143e-01, -2.21844482e-01,  1.84997223e+00,
      ...
              -1.21051799e+00, -3.67753776e-01,  1.00004975e+00,
              -1.30402450e-01,  1.08847436e-01, -8.06239420e-01,
               3.04261341e-01, -2.75164749e-01, -5.39020456e-01,
              -1.07610620e-01,  2.18530560e+00,  4.41874710e-01,
              -2.37728559e-01, -4.50805654e-01,  2.72418581e-01,
              -7.38092586e-01, -2.34393681e+00,  5.21358491e-01,
              -1.02342427e+00, -2.27437571e-01,  8.34447225e-01,
               5.25093339e-01, -9.74835311e-01,  1.28344006e+00,
               1.31464981e+00, -1.21745550e-01, -2.18126345e-01,
               2.30982110e-01, -1.40217342e+00,  6.51477656e-01,
               8.27257637e-01, -1.74879594e-01,  1.34584346e+00,
               1.35776443e+00, -1.46789158e+00,  9.23491710e-02,
              -1.54262898e-01,  9.00426851e-01, -9.76272546e-01,
               5.84051037e-01, -7.20015378e-02,  3.90502365e-01,
              -1.00072682e+00,  6.26765863e-01, -6.57189342e-01,
               5.71041199e-01, -4.38164033e-01, -1.44436373e+00,
               4.67760783e-01,  1.38629530e+00,  3.73237264e-01,
              -9.37224147e-01,  1.95161086e-01, -9.65195629e-01,
              -1.61252987e+00,  3.06965268e-01, -1.25784537e-01,
              -1.08614370e+00, -3.63936865e-02]])

    • p
      (chain, draw, obs, alts_probs)
      float64
      1.347e-100 1.911e-14 ... 1.0 0.0

      array([[[[1.34671297e-100, 1.91139328e-014, 2.43947490e-063,
                1.00000000e+000, 0.00000000e+000],
               [9.98040992e-063, 1.00000000e+000, 1.02487442e-070,
                4.07132826e-052, 0.00000000e+000],
               [2.85515600e-112, 1.00000000e+000, 1.30621444e-158,
                1.08048464e-045, 0.00000000e+000],
               ...,
               [2.30156887e-160, 3.08592709e-038, 5.83005298e-165,
                1.00000000e+000, 0.00000000e+000],
               [6.06929797e-204, 9.61695016e-001, 1.99988451e-140,
                3.83049845e-002, 0.00000000e+000],
               [3.26382940e-124, 4.08947632e-044, 4.62453374e-136,
                1.00000000e+000, 0.00000000e+000]],
      
              [[1.00000000e+000, 1.41883359e-028, 1.25084719e-174,
                2.54756743e-201, 1.13562547e-244],
               [1.00000000e+000, 1.79329481e-020, 3.48460260e-171,
                3.46461278e-172, 1.86502610e-230],
               [1.00000000e+000, 2.02385296e-023, 3.71931404e-131,
                9.58177773e-151, 4.49320481e-193],
      ...
               [0.00000000e+000, 0.00000000e+000, 0.00000000e+000,
                0.00000000e+000, 1.00000000e+000],
               [0.00000000e+000, 0.00000000e+000, 0.00000000e+000,
                0.00000000e+000, 1.00000000e+000],
               [0.00000000e+000, 0.00000000e+000, 0.00000000e+000,
                0.00000000e+000, 1.00000000e+000]],
      
              [[1.98431291e-105, 1.00000000e+000, 3.58445356e-095,
                1.78993656e-020, 0.00000000e+000],
               [3.84818059e-076, 1.00000000e+000, 2.71493217e-123,
                1.32300540e-099, 0.00000000e+000],
               [3.51959229e-139, 1.00000000e+000, 5.82198856e-227,
                4.52246690e-086, 0.00000000e+000],
               ...,
               [1.42197003e-168, 3.25818873e-016, 2.21723790e-206,
                1.00000000e+000, 0.00000000e+000],
               [1.12932342e-249, 1.00000000e+000, 1.15561543e-198,
                7.89297474e-029, 0.00000000e+000],
               [2.17119127e-122, 1.05514439e-026, 1.52717246e-173,
                1.00000000e+000, 0.00000000e+000]]]])

  • Indexes: (4)
    • chain
      PandasIndex

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

    • draw
      PandasIndex

      PandasIndex(Index([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,
             ...
             490, 491, 492, 493, 494, 495, 496, 497, 498, 499],
            dtype='int64', name='draw', length=500))

    • obs
      PandasIndex

      PandasIndex(Index([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,
             ...
             890, 891, 892, 893, 894, 895, 896, 897, 898, 899],
            dtype='int64', name='obs', length=900))

    • alts_probs
      PandasIndex

      PandasIndex(Index(['ec', 'er', 'gc', 'gr', 'hp'], dtype='object', name='alts_probs'))

  • Attributes: (4)

    created_at :

    2024-06-23T20:32:10.756506+00:00

    arviz_version :

    0.18.0

    inference_library :

    pymc

    inference_library_version :

    5.15.1

  
  
• prior_predictive

  <xarray.Dataset> Size: 4MB
  Dimensions:  (chain: 1, draw: 500, obs: 900)
  Coordinates:
    * chain    (chain) int64 8B 0
    * draw     (draw) int64 4kB 0 1 2 3 4 5 6 7 ... 493 494 495 496 497 498 499
    * obs      (obs) int64 7kB 0 1 2 3 4 5 6 7 ... 892 893 894 895 896 897 898 899
  Data variables:
      y_cat    (chain, draw, obs) int64 4MB 3 1 1 2 1 3 3 3 3 ... 1 1 1 3 1 3 1 3
  Attributes:
      created_at:                 2024-06-23T20:32:10.758790+00:00
      arviz_version:              0.18.0
      inference_library:          pymc
      inference_library_version:  5.15.1

  xarray.Dataset

  • Dimensions:
    • chain: 1
    • draw: 500
    • obs: 900
  • Coordinates: (3)
    • chain
      (chain)
      int64
      0

      array([0])

    • draw
      (draw)
      int64
      0 1 2 3 4 5 ... 495 496 497 498 499

      array([  0,   1,   2, ..., 497, 498, 499])

    • obs
      (obs)
      int64
      0 1 2 3 4 5 ... 895 896 897 898 899

      array([  0,   1,   2, ..., 897, 898, 899])

  • Data variables: (1)
    • y_cat
      (chain, draw, obs)
      int64
      3 1 1 2 1 3 3 3 ... 1 1 1 3 1 3 1 3

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

  • Indexes: (3)
    • chain
      PandasIndex

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

    • draw
      PandasIndex

      PandasIndex(Index([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,
             ...
             490, 491, 492, 493, 494, 495, 496, 497, 498, 499],
            dtype='int64', name='draw', length=500))

    • obs
      PandasIndex

      PandasIndex(Index([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,
             ...
             890, 891, 892, 893, 894, 895, 896, 897, 898, 899],
            dtype='int64', name='obs', length=900))

  • Attributes: (4)

    created_at :

    2024-06-23T20:32:10.758790+00:00

    arviz_version :

    0.18.0

    inference_library :

    pymc

    inference_library_version :

    5.15.1

  
  
• observed_data

  <xarray.Dataset> Size: 14kB
  Dimensions:  (obs: 900)
  Coordinates:
    * obs      (obs) int64 7kB 0 1 2 3 4 5 6 7 ... 892 893 894 895 896 897 898 899
  Data variables:
      y_cat    (obs) int64 7kB 2 2 2 1 1 2 2 2 2 2 2 2 ... 2 2 2 2 2 2 1 2 2 2 2 2
  Attributes:
      created_at:                 2024-06-23T20:32:10.759288+00:00
      arviz_version:              0.18.0
      inference_library:          pymc
      inference_library_version:  5.15.1

  xarray.Dataset

  • Dimensions:
    • obs: 900
  • Coordinates: (1)
    • obs
      (obs)
      int64
      0 1 2 3 4 5 ... 895 896 897 898 899

      array([  0,   1,   2, ..., 897, 898, 899])

  • Data variables: (1)
    • y_cat
      (obs)
      int64
      2 2 2 1 1 2 2 2 ... 2 2 1 2 2 2 2 2

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

  • Indexes: (1)
    • obs
      PandasIndex

      PandasIndex(Index([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,
             ...
             890, 891, 892, 893, 894, 895, 896, 897, 898, 899],
            dtype='int64', name='obs', length=900))

  • Attributes: (4)

    created_at :

    2024-06-23T20:32:10.759288+00:00

    arviz_version :

    0.18.0

    inference_library :

    pymc

    inference_library_version :

    5.15.1

  
  

In [7]:

summaries = az.summary(idata_m1, var_names=["beta_ic", "beta_oc"])
summaries # results are sus because beta_ic and beta_oc should both be negative

	mean	sd	hdi_3%	hdi_97%	mcse_mean	mcse_sd	ess_bulk	ess_tail	r_hat
beta_ic	0.002	0.0	0.002	0.002	0.0	0.0	1802.0	2209.0	1.0
beta_oc	-0.004	0.0	-0.005	-0.004	0.0	0.0	983.0	1241.0	1.0

In [8]:

## marginal rate of substitution for a reduction in installation costs
post = az.extract(idata_m1)
substitution_rate = post["beta_oc"] / post["beta_ic"]
substitution_rate.mean().item()

-2.211196691253724

In [9]:

fig, ax = plt.subplots(figsize=(20, 10))

ax.hist(
    substitution_rate,
    bins=30,
    ec="black",
)
ax.set_title("Uncertainty in Marginal Rate of Substitution \n Operating Costs / Installation Costs");
plt.show()

In [10]:

idata_m1["posterior"]["p"].mean(dim=["chain", "draw", "obs"])

<xarray.DataArray 'p' (alts_probs: 5)> Size: 40B
array([0.08432415, 0.13772203, 0.26912233, 0.38199249, 0.12683901])
Coordinates:
  * alts_probs  (alts_probs) <U2 40B 'ec' 'er' 'gc' 'gr' 'hp'

xarray.DataArray

'p'

• alts_probs: 5

• 0.08432 0.1377 0.2691 0.382 0.1268

  array([0.08432415, 0.13772203, 0.26912233, 0.38199249, 0.12683901])

• Coordinates: (1)
  • alts_probs
    (alts_probs)
    <U2
    'ec' 'er' 'gc' 'gr' 'hp'

    array(['ec', 'er', 'gc', 'gr', 'hp'], dtype='<U2')

• Indexes: (1)
  • alts_probs
    PandasIndex

    PandasIndex(Index(['ec', 'er', 'gc', 'gr', 'hp'], dtype='object', name='alts_probs'))

• Attributes: (0)

In [11]:

fig, axs = plt.subplots(1, 2, figsize=(20, 10))
ax = axs[0]
counts = wide_heating_df.groupby("depvar")["idcase"].count()
predicted_shares = idata_m1["posterior"]["p"].mean(dim=["chain", "draw", "obs"])
ci_lb = idata_m1["posterior"]["p"].quantile(0.025, dim=["chain", "draw", "obs"])
ci_ub = idata_m1["posterior"]["p"].quantile(0.975, dim=["chain", "draw", "obs"])
ax.scatter(ci_lb, ["ec", "er", "gc", "gr", "hp"], color="k", s=2)
ax.scatter(ci_ub, ["ec", "er", "gc", "gr", "hp"], color="k", s=2)
ax.scatter(
    counts / counts.sum(),
    ["ec", "er", "gc", "gr", "hp"],
    label="Observed Shares",
    color="red",
    s=100,
)
ax.hlines(
    ["ec", "er", "gc", "gr", "hp"], ci_lb, ci_ub, label="Predicted 95% Interval", color="black"
)
ax.legend()
ax.set_title("Observed V Predicted Shares")
az.plot_ppc(idata_m1, ax=axs[1])
axs[1].set_title("Posterior Predictive Checks")
ax.set_xlabel("Shares")
ax.set_ylabel("Heating System");

Getting a better fit

You can improve the fit by adding a good-specific intercept – intuitively this should allow you to match the observed frequencies

In [12]:

N = wide_heating_df.shape[0]
observed = pd.Categorical(wide_heating_df["depvar"]).codes

coords = {
    "alts_intercepts": ["ec", "er", "gc", "gr"],
    "alts_probs": ["ec", "er", "gc", "gr", "hp"],
    "obs": range(N),
}
with pm.Model(coords=coords) as model_2:
    beta_ic = pm.Normal("beta_ic", 0, 1)
    beta_oc = pm.Normal("beta_oc", 0, 1)
    alphas = pm.Normal("alpha", 0, 1, dims="alts_intercepts")

    ## Construct Utility matrix and Pivot using an intercept per alternative
    u0 = alphas[0] + beta_ic * wide_heating_df["ic.ec"] + beta_oc * wide_heating_df["oc.ec"]
    u1 = alphas[1] + beta_ic * wide_heating_df["ic.er"] + beta_oc * wide_heating_df["oc.er"]
    u2 = alphas[2] + beta_ic * wide_heating_df["ic.gc"] + beta_oc * wide_heating_df["oc.gc"]
    u3 = alphas[3] + beta_ic * wide_heating_df["ic.gr"] + beta_oc * wide_heating_df["oc.gr"]
    u4 = np.zeros(N)  # Outside Good
    s = pm.math.stack([u0, u1, u2, u3, u4]).T

    ## Apply Softmax Transform
    p_ = pm.Deterministic("p", pm.math.softmax(s, axis=1), dims=("obs", "alts_probs"))

    ## Likelihood
    choice_obs = pm.Categorical("y_cat", p=p_, observed=observed, dims="obs")

    idata_m2 = pm.sample_prior_predictive()
    idata_m2.extend(
        pm.sample(nuts_sampler="numpyro", idata_kwargs={"log_likelihood": True}, random_seed=103)
    )
    idata_m2.extend(pm.sample_posterior_predictive(idata_m2))


pm.model_to_graphviz(model_2)

Sampling: [alpha, beta_ic, beta_oc, y_cat]

Sampling: [y_cat]

In [13]:

az.summary(idata_m2, var_names=["beta_ic", "beta_oc", "alpha"])

# model is still terrible, installation costs and operational costs are meaningless, it's all about the intercepts

	mean	sd	hdi_3%	hdi_97%	mcse_mean	mcse_sd	ess_bulk	ess_tail	r_hat
beta_ic	0.001	0.000	-0.000	0.002	0.000	0.000	1324.0	1904.0	1.00
beta_oc	-0.003	0.001	-0.005	-0.001	0.000	0.000	1337.0	1785.0	1.00
alpha[ec]	1.063	0.513	0.087	2.054	0.017	0.012	914.0	1150.0	1.00
alpha[er]	1.102	0.491	0.135	2.017	0.017	0.012	856.0	900.0	1.01
alpha[gc]	2.397	0.320	1.826	3.019	0.011	0.008	874.0	1002.0	1.01
alpha[gr]	0.760	0.385	0.036	1.448	0.012	0.009	957.0	1124.0	1.01

In [14]:

fig, axs = plt.subplots(1, 2, figsize=(20, 10))
ax = axs[0]
counts = wide_heating_df.groupby("depvar")["idcase"].count()
predicted_shares = idata_m2["posterior"]["p"].mean(dim=["chain", "draw", "obs"])
ci_lb = idata_m2["posterior"]["p"].quantile(0.025, dim=["chain", "draw", "obs"])
ci_ub = idata_m2["posterior"]["p"].quantile(0.975, dim=["chain", "draw", "obs"])

ax.scatter(ci_lb, ["ec", "er", "gc", "gr", "hp"], color="k", s=2)
ax.scatter(ci_ub, ["ec", "er", "gc", "gr", "hp"], color="k", s=2)
ax.scatter(
    counts / counts.sum(),
    ["ec", "er", "gc", "gr", "hp"],
    label="Observed Shares",
    color="red",
    s=100,
)
ax.hlines(
    ["ec", "er", "gc", "gr", "hp"], ci_lb, ci_ub, label="Predicted 95% Interval", color="black"
)
ax.legend()
ax.set_title("Observed V Predicted Shares")
az.plot_ppc(idata_m2, ax=axs[1])
axs[1].set_title("Posterior Predictive Checks")
ax.set_xlabel("Shares")
ax.set_ylabel("Heating System");

In [15]:

## marginal rate of substitution for a reduction in installation costs
post = az.extract(idata_m2)
substitution_rate = post["beta_oc"] / post["beta_ic"]
substitution_rate.mean().item()

-11.215830921463027

In [16]:

fig, ax = plt.subplots(figsize=(20, 10))

ax.hist(
    substitution_rate,
    bins=30,
    ec="black",
)
ax.set_title("Uncertainty in Marginal Rate of Substitution \n Operating Costs / Installation Costs");

Bells and Whistles

Adding some cool extra features by: - using a multivariate normal distirbution – this allows for correlation among alternatives - controlling for income in a log-affine fashion

Intuitively not sure how the correlation is identified, given that you only observe the actual choice

In [17]:

coords = {
    "alts_intercepts": ["ec", "er", "gc", "gr"],
    "alts_probs": ["ec", "er", "gc", "gr", "hp"],
    "obs": range(N),
}
with pm.Model(coords=coords) as model_3:
    ## Add data to experiment with changes later.
    ic_ec = pm.Data("ic_ec", wide_heating_df["ic.ec"])
    oc_ec = pm.Data("oc_ec", wide_heating_df["oc.ec"])
    ic_er = pm.Data("ic_er", wide_heating_df["ic.er"])
    oc_er = pm.Data("oc_er", wide_heating_df["oc.er"])

    beta_ic = pm.Normal("beta_ic", 0, 1)
    beta_oc = pm.Normal("beta_oc", 0, 1)
    beta_income = pm.Normal("beta_income", 0, 1, dims="alts_intercepts")
    chol, corr, stds = pm.LKJCholeskyCov(
        "chol", n=4, eta=2.0, sd_dist=pm.Exponential.dist(1.0, shape=4)
    )
    alphas = pm.MvNormal("alpha", mu=0, chol=chol, dims="alts_intercepts")

    u0 = alphas[0] + beta_ic * ic_ec + beta_oc * oc_ec + beta_income[0] * wide_heating_df["income"]
    u1 = alphas[1] + beta_ic * ic_er + beta_oc * oc_er + beta_income[1] * wide_heating_df["income"]
    u2 = (
        alphas[2]
        + beta_ic * wide_heating_df["ic.gc"]
        + beta_oc * wide_heating_df["oc.gc"]
        + beta_income[2] * wide_heating_df["income"]
    )
    u3 = (
        alphas[3]
        + beta_ic * wide_heating_df["ic.gr"]
        + beta_oc * wide_heating_df["oc.gr"]
        + beta_income[3] * wide_heating_df["income"]
    )
    u4 = np.zeros(N)  # pivot
    s = pm.math.stack([u0, u1, u2, u3, u4]).T

    p_ = pm.Deterministic("p", pm.math.softmax(s, axis=1), dims=("obs", "alts_probs"))
    choice_obs = pm.Categorical("y_cat", p=p_, observed=observed, dims="obs")

    idata_m3 = pm.sample_prior_predictive()
    idata_m3.extend(
        pm.sample(nuts_sampler="numpyro", idata_kwargs={"log_likelihood": True}, random_seed=100)
    )
    idata_m3.extend(pm.sample_posterior_predictive(idata_m3))


pm.model_to_graphviz(model_3)

Sampling: [alpha, beta_ic, beta_income, beta_oc, chol, y_cat]

There were 713 divergences after tuning. Increase `target_accept` or reparameterize.
The rhat statistic is larger than 1.01 for some parameters. This indicates problems during sampling. See https://arxiv.org/abs/1903.08008 for details
The effective sample size per chain is smaller than 100 for some parameters.  A higher number is needed for reliable rhat and ess computation. See https://arxiv.org/abs/1903.08008 for details
Sampling: [y_cat]

In [18]:

fig, axs = plt.subplots(1, 2, figsize=(20, 10))
ax = axs[0]
counts = wide_heating_df.groupby("depvar")["idcase"].count()
predicted_shares = idata_m3["posterior"]["p"].mean(dim=["chain", "draw", "obs"])
ci_lb = idata_m3["posterior"]["p"].quantile(0.025, dim=["chain", "draw", "obs"])
ci_ub = idata_m3["posterior"]["p"].quantile(0.975, dim=["chain", "draw", "obs"])

ax.scatter(ci_lb, ["ec", "er", "gc", "gr", "hp"], color="k", s=2)
ax.scatter(ci_ub, ["ec", "er", "gc", "gr", "hp"], color="k", s=2)
ax.scatter(
    counts / counts.sum(),
    ["ec", "er", "gc", "gr", "hp"],
    label="Observed Shares",
    color="red",
    s=100,
)
ax.hlines(
    ["ec", "er", "gc", "gr", "hp"], ci_lb, ci_ub, label="Predicted 95% Interval", color="black"
)
ax.legend()
ax.set_title("Observed V Predicted Shares")
az.plot_ppc(idata_m3, ax=axs[1])
axs[1].set_title("Posterior Predictive Checks")
ax.set_xlabel("Shares")
ax.set_ylabel("Heating System");

In [19]:

az.summary(
    idata_m3, var_names=["beta_income", "beta_ic", "beta_oc", "alpha", "chol_corr"], round_to=4
)

/Users/gzheng/projects/professional/likert_dualresponse/.venv/lib/python3.11/site-packages/arviz/stats/diagnostics.py:592: RuntimeWarning: invalid value encountered in scalar divide
  (between_chain_variance / within_chain_variance + num_samples - 1) / (num_samples)

	mean	sd	hdi_3%	hdi_97%	mcse_mean	mcse_sd	ess_bulk	ess_tail	r_hat
beta_income[ec]	0.1090	0.1045	-0.1194	0.2745	0.0111	0.0079	125.5245	1785.3613	1.0326
beta_income[er]	0.0568	0.0985	-0.1462	0.2349	0.0034	0.0081	814.4173	1351.8588	1.1603
beta_income[gc]	0.0701	0.0810	-0.0908	0.2230	0.0039	0.0028	514.8539	1447.0504	1.0928
beta_income[gr]	-0.0182	0.0991	-0.2034	0.1530	0.0166	0.0119	43.2790	1394.2747	1.0660
beta_ic	0.0005	0.0007	-0.0008	0.0016	0.0002	0.0001	21.2834	1025.6984	1.1204
beta_oc	-0.0037	0.0014	-0.0063	-0.0008	0.0001	0.0001	130.7170	1819.1965	1.0273
alpha[ec]	0.9330	1.0135	-0.4580	2.9612	0.1099	0.0780	361.4979	1033.9501	1.0793
alpha[er]	1.1878	1.0135	-0.3409	3.1845	0.0761	0.0608	377.7638	1058.1381	1.1560
alpha[gc]	2.2653	0.7507	1.2054	3.7209	0.1258	0.0897	46.5848	959.9773	1.0614
alpha[gr]	1.0156	0.9110	-0.1379	2.6814	0.2079	0.1493	17.7565	40.8513	1.1541
chol_corr[0, 0]	1.0000	0.0000	1.0000	1.0000	0.0000	0.0000	4000.0000	4000.0000	NaN
chol_corr[0, 1]	0.1392	0.3365	-0.5441	0.7382	0.0075	0.0220	2135.5951	2551.8261	1.1486
chol_corr[0, 2]	0.0750	0.3488	-0.5012	0.7344	0.0588	0.0419	47.6418	1379.3299	1.0662
chol_corr[0, 3]	0.0536	0.3764	-0.5126	0.7631	0.0860	0.0618	24.2892	2253.6964	1.1096
chol_corr[1, 0]	0.1392	0.3365	-0.5441	0.7382	0.0075	0.0220	2135.5951	2551.8261	1.1486
chol_corr[1, 1]	1.0000	0.0000	1.0000	1.0000	0.0000	0.0000	38.4914	2944.7704	1.0696
chol_corr[1, 2]	0.1438	0.3286	-0.4194	0.8010	0.0223	0.0344	199.9289	2081.5502	1.0314
chol_corr[1, 3]	0.0821	0.3555	-0.4699	0.7795	0.0605	0.0431	45.2660	2019.6281	1.0575
chol_corr[2, 0]	0.0750	0.3488	-0.5012	0.7344	0.0588	0.0419	47.6418	1379.3299	1.0662
chol_corr[2, 1]	0.1438	0.3286	-0.4194	0.8010	0.0223	0.0344	199.9289	2081.5502	1.0314
chol_corr[2, 2]	1.0000	0.0000	1.0000	1.0000	0.0000	0.0000	21.0671	3321.0770	1.1213
chol_corr[2, 3]	0.0885	0.4025	-0.5738	0.6999	0.1268	0.0924	12.2708	4.3794	1.2387
chol_corr[3, 0]	0.0536	0.3764	-0.5126	0.7631	0.0860	0.0618	24.2892	2253.6964	1.1096
chol_corr[3, 1]	0.0821	0.3555	-0.4699	0.7795	0.0605	0.0431	45.2660	2019.6281	1.0575
chol_corr[3, 2]	0.0885	0.4025	-0.5738	0.6999	0.1268	0.0924	12.2708	4.3794	1.2387
chol_corr[3, 3]	1.0000	0.0000	1.0000	1.0000	0.0000	0.0000	3464.6160	3626.8755	1.0597

In [20]:

post = az.extract(idata_m3)
substitution_rate = post["beta_oc"] / post["beta_ic"]
substitution_rate.mean().item()

-13.242388002273604

Tweaks

My preferred model is model 2, but with - a control for income - (ideally) set-constrained coefficients for installation cost and operating cost. However, this doesn’t work. In practice it’s because the coefficients are both so small (especially the \(\beta_\text{installation cost}\)) that constraining it is difficult

Because they’re so small and so close to zero, it’s also not meaningful to talk about a MRS.

In [21]:

coords = {
    "alts_intercepts": ["ec", "er", "gc", "gr"],
    "alts_probs": ["ec", "er", "gc", "gr", "hp"],
    "obs": range(N),
}

with pm.Model(coords=coords) as model_4:
    beta_ic = pm.Normal("beta_ic")
    beta_oc = pm.Normal("beta_oc")
    beta_income = pm.Normal("beta_income", 0, 1, dims="alts_intercepts")
    alphas = pm.Normal("alpha", 0, 5, dims="alts_intercepts")

    ## Construct Utility matrix and Pivot using an intercept per alternative
    u0 = alphas[0] + beta_ic * wide_heating_df["ic.ec"] + beta_oc * wide_heating_df["oc.ec"] + beta_income[0] * wide_heating_df["income"]
    u1 = alphas[1] + beta_ic * wide_heating_df["ic.er"] + beta_oc * wide_heating_df["oc.er"] + beta_income[1] * wide_heating_df["income"]
    u2 = alphas[2] + beta_ic * wide_heating_df["ic.gc"] + beta_oc * wide_heating_df["oc.gc"] + beta_income[2] * wide_heating_df["income"]
    u3 = alphas[3] + beta_ic * wide_heating_df["ic.gr"] + beta_oc * wide_heating_df["oc.gr"] + beta_income[3] * wide_heating_df["income"]
    u4 = np.zeros(N)  # Outside Good
    s = pm.math.stack([u0, u1, u2, u3, u4]).T

    ## Apply Softmax Transform
    p_ = pm.Deterministic("p", pm.math.softmax(s, axis=1), dims=("obs", "alts_probs"))

    ## Likelihood
    choice_obs = pm.Categorical("y_cat", p=p_, observed=observed, dims="obs")

    idata_m4 = pm.sample_prior_predictive()
    idata_m4.extend(
        pm.sample(nuts_sampler="numpyro", idata_kwargs={"log_likelihood": True}, random_seed=100)
    )
    idata_m4.extend(pm.sample_posterior_predictive(idata_m4))


pm.model_to_graphviz(model_4)

Sampling: [alpha, beta_ic, beta_income, beta_oc, y_cat]

Sampling: [y_cat]

In [22]:

fig, axs = plt.subplots(1, 2, figsize=(20, 10))
ax = axs[0]
counts = wide_heating_df.groupby("depvar")["idcase"].count()
predicted_shares = idata_m4["posterior"]["p"].mean(dim=["chain", "draw", "obs"])
ci_lb = idata_m4["posterior"]["p"].quantile(0.025, dim=["chain", "draw", "obs"])
ci_ub = idata_m4["posterior"]["p"].quantile(0.975, dim=["chain", "draw", "obs"])

ax.scatter(ci_lb, ["ec", "er", "gc", "gr", "hp"], color="k", s=2)
ax.scatter(ci_ub, ["ec", "er", "gc", "gr", "hp"], color="k", s=2)
ax.scatter(
    counts / counts.sum(),
    ["ec", "er", "gc", "gr", "hp"],
    label="Observed Shares",
    color="red",
    s=100,
)
ax.hlines(
    ["ec", "er", "gc", "gr", "hp"], ci_lb, ci_ub, label="Predicted 95% Interval", color="black"
)
ax.legend()
ax.set_title("Observed V Predicted Shares")
az.plot_ppc(idata_m4, ax=axs[1])
axs[1].set_title("Posterior Predictive Checks")
ax.set_xlabel("Shares")
ax.set_ylabel("Heating System");

In [23]:

az.summary(
    idata_m4, var_names=["beta_ic", "beta_oc", "beta_income", "alpha"], round_to=4
)

	mean	sd	hdi_3%	hdi_97%	mcse_mean	mcse_sd	ess_bulk	ess_tail	r_hat
beta_ic	-0.0007	0.0006	-0.0018	0.0004	0.0000	0.0000	2034.8098	2619.1608	1.0010
beta_oc	-0.0055	0.0015	-0.0081	-0.0025	0.0000	0.0000	2549.3952	3052.4259	1.0023
beta_income[ec]	-0.0415	0.1123	-0.2552	0.1597	0.0032	0.0023	1233.5860	1949.2437	1.0024
beta_income[er]	-0.0708	0.1064	-0.2544	0.1412	0.0031	0.0022	1151.0759	1873.4452	1.0022
beta_income[gc]	-0.0527	0.0866	-0.2195	0.1042	0.0027	0.0019	1063.0967	1595.1725	1.0012
beta_income[gr]	-0.1619	0.0990	-0.3460	0.0231	0.0029	0.0021	1134.4858	1809.5429	1.0011
alpha[ec]	3.5846	0.9338	1.9701	5.4876	0.0286	0.0204	1066.7882	1428.9757	1.0018
alpha[er]	3.8573	0.9128	2.1368	5.5381	0.0294	0.0210	968.1931	1275.9228	1.0013
alpha[gc]	4.1892	0.6397	2.9483	5.3714	0.0211	0.0150	918.4178	1199.7850	1.0012
alpha[gr]	3.1823	0.7224	1.8973	4.5941	0.0232	0.0165	967.8811	1118.1689	1.0017

In [24]:

post = az.extract(idata_m4)
substitution_rate = post["beta_oc"] / post["beta_ic"]
substitution_rate.mean().item()

31.17482352077738

In [25]:

fig, ax = plt.subplots(figsize=(20, 10))

ax.hist(
    substitution_rate,
    bins=30,
    ec="black",
)
ax.set_title("Uncertainty in Marginal Rate of Substitution \n Operating Costs / Installation Costs");

Counterfactuals

In [26]:

with model_3:
    # update values of predictors with new 20% price increase in operating costs for electrical options
    pm.set_data({"oc_ec": wide_heating_df["oc.ec"] * 1.2, "oc_er": wide_heating_df["oc.er"] * 1.2})
    # use the updated values and predict outcomes and probabilities:
    idata_new_policy = pm.sample_posterior_predictive(
        idata_m3,
        var_names=["p", "y_cat"],
        return_inferencedata=True,
        predictions=True,
        extend_inferencedata=False,
        random_seed=100,
    )

idata_new_policy

Sampling: [y_cat]

arviz.InferenceData

• predictions

  <xarray.Dataset> Size: 173MB
  Dimensions:     (chain: 4, draw: 1000, obs: 900, alts_probs: 5)
  Coordinates:
    * chain       (chain) int64 32B 0 1 2 3
    * draw        (draw) int64 8kB 0 1 2 3 4 5 6 7 ... 993 994 995 996 997 998 999
    * obs         (obs) int64 7kB 0 1 2 3 4 5 6 7 ... 893 894 895 896 897 898 899
    * alts_probs  (alts_probs) <U2 40B 'ec' 'er' 'gc' 'gr' 'hp'
  Data variables:
      p           (chain, draw, obs, alts_probs) float64 144MB 0.05725 ... 0.04717
      y_cat       (chain, draw, obs) int64 29MB 2 4 4 2 1 2 2 2 ... 3 3 2 1 2 4 2
  Attributes:
      created_at:                 2024-06-23T20:34:06.401295+00:00
      arviz_version:              0.18.0
      inference_library:          pymc
      inference_library_version:  5.15.1

  xarray.Dataset

  • Dimensions:
    • chain: 4
    • draw: 1000
    • obs: 900
    • alts_probs: 5
  • Coordinates: (4)
    • 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
      (obs)
      int64
      0 1 2 3 4 5 ... 895 896 897 898 899

      array([  0,   1,   2, ..., 897, 898, 899])

    • alts_probs
      (alts_probs)
      <U2
      'ec' 'er' 'gc' 'gr' 'hp'

      array(['ec', 'er', 'gc', 'gr', 'hp'], dtype='<U2')

  • Data variables: (2)
    • p
      (chain, draw, obs, alts_probs)
      float64
      0.05725 0.06145 ... 0.1894 0.04717

      array([[[[0.05724768, 0.06144754, 0.68573654, 0.12073992, 0.07482832],
               [0.0536941 , 0.06504792, 0.66728469, 0.13613161, 0.07784168],
               [0.06074645, 0.07880667, 0.63082648, 0.15196209, 0.0776583 ],
               ...,
               [0.05956621, 0.07063259, 0.65266042, 0.13936663, 0.07777414],
               [0.06049453, 0.0836567 , 0.63028224, 0.14349363, 0.0820729 ],
               [0.04665807, 0.07612694, 0.61031866, 0.17417927, 0.09271705]],
      
              [[0.0575523 , 0.05957472, 0.70580507, 0.12394815, 0.05311976],
               [0.04978846, 0.05194444, 0.70624589, 0.13173661, 0.0602846 ],
               [0.05858354, 0.06523037, 0.65436287, 0.15864973, 0.06317349],
               ...,
               [0.0572611 , 0.06067699, 0.67328214, 0.15106374, 0.05771602],
               [0.05613799, 0.08006781, 0.64611688, 0.15639254, 0.06128477],
               [0.03569163, 0.04690307, 0.65923192, 0.18508592, 0.07308747]],
      
              [[0.03848548, 0.06967034, 0.71777804, 0.11688472, 0.05718142],
               [0.04106635, 0.06133217, 0.70216202, 0.12892532, 0.06651414],
               [0.05499752, 0.08312547, 0.62147335, 0.17002871, 0.07037495],
               ...,
      ...
               ...,
               [0.05812185, 0.07607232, 0.66578602, 0.14396336, 0.05605645],
               [0.06088941, 0.08638921, 0.64408592, 0.14895731, 0.05967816],
               [0.05234028, 0.06589006, 0.65931254, 0.14420871, 0.07824841]],
      
              [[0.04946695, 0.04037203, 0.74481053, 0.10060332, 0.06474716],
               [0.05557851, 0.04546533, 0.71803102, 0.10593543, 0.07498971],
               [0.07094328, 0.07473905, 0.60520828, 0.16683318, 0.08227621],
               ...,
               [0.05864539, 0.05887554, 0.64364577, 0.17843161, 0.0604017 ],
               [0.04773015, 0.08160192, 0.63828546, 0.17464675, 0.05773572],
               [0.04628499, 0.05529322, 0.60885068, 0.24491774, 0.04465337]],
      
              [[0.03588243, 0.06452423, 0.70246049, 0.13061053, 0.06652232],
               [0.04836715, 0.064673  , 0.69243603, 0.12374824, 0.07077557],
               [0.0621403 , 0.08466819, 0.62086437, 0.15856235, 0.07376478],
               ...,
               [0.04919897, 0.07558963, 0.6408676 , 0.17349437, 0.06084943],
               [0.04264316, 0.09181579, 0.6372937 , 0.16979284, 0.05845452],
               [0.06034545, 0.06557241, 0.6375511 , 0.18935961, 0.04717142]]]])

    • y_cat
      (chain, draw, obs)
      int64
      2 4 4 2 1 2 2 2 ... 3 3 3 2 1 2 4 2

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

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

      PandasIndex(Index([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,
             ...
             890, 891, 892, 893, 894, 895, 896, 897, 898, 899],
            dtype='int64', name='obs', length=900))

    • alts_probs
      PandasIndex

      PandasIndex(Index(['ec', 'er', 'gc', 'gr', 'hp'], dtype='object', name='alts_probs'))

  • Attributes: (4)

    created_at :

    2024-06-23T20:34:06.401295+00:00

    arviz_version :

    0.18.0

    inference_library :

    pymc

    inference_library_version :

    5.15.1

  
  
• predictions_constant_data

  <xarray.Dataset> Size: 58kB
  Dimensions:      (ic_ec_dim_0: 900, ic_er_dim_0: 900, oc_ec_dim_0: 900,
                    oc_er_dim_0: 900)
  Coordinates:
    * ic_ec_dim_0  (ic_ec_dim_0) int64 7kB 0 1 2 3 4 5 ... 894 895 896 897 898 899
    * ic_er_dim_0  (ic_er_dim_0) int64 7kB 0 1 2 3 4 5 ... 894 895 896 897 898 899
    * oc_ec_dim_0  (oc_ec_dim_0) int64 7kB 0 1 2 3 4 5 ... 894 895 896 897 898 899
    * oc_er_dim_0  (oc_er_dim_0) int64 7kB 0 1 2 3 4 5 ... 894 895 896 897 898 899
  Data variables:
      ic_ec        (ic_ec_dim_0) float64 7kB 859.9 796.8 719.9 ... 799.8 967.9
      ic_er        (ic_er_dim_0) float64 7kB 995.8 894.7 ... 1.123e+03 1.092e+03
      oc_ec        (oc_ec_dim_0) float64 7kB 664.0 624.3 526.9 ... 594.2 622.4
      oc_er        (oc_er_dim_0) float64 7kB 606.7 583.8 485.7 ... 481.9 550.2
  Attributes:
      created_at:                 2024-06-23T20:34:06.404039+00:00
      arviz_version:              0.18.0
      inference_library:          pymc
      inference_library_version:  5.15.1

  xarray.Dataset

  • Dimensions:
    • ic_ec_dim_0: 900
    • ic_er_dim_0: 900
    • oc_ec_dim_0: 900
    • oc_er_dim_0: 900
  • Coordinates: (4)
    • ic_ec_dim_0
      (ic_ec_dim_0)
      int64
      0 1 2 3 4 5 ... 895 896 897 898 899

      array([  0,   1,   2, ..., 897, 898, 899])

    • ic_er_dim_0
      (ic_er_dim_0)
      int64
      0 1 2 3 4 5 ... 895 896 897 898 899

      array([  0,   1,   2, ..., 897, 898, 899])

    • oc_ec_dim_0
      (oc_ec_dim_0)
      int64
      0 1 2 3 4 5 ... 895 896 897 898 899

      array([  0,   1,   2, ..., 897, 898, 899])

    • oc_er_dim_0
      (oc_er_dim_0)
      int64
      0 1 2 3 4 5 ... 895 896 897 898 899

      array([  0,   1,   2, ..., 897, 898, 899])

  • Data variables: (4)
    • ic_ec
      (ic_ec_dim_0)
      float64
      859.9 796.8 719.9 ... 799.8 967.9

      array([ 859.9 ,  796.82,  719.86,  761.25,  858.86,  693.74,  633.63,
              813.33,  875.7 ,  775.86,  792.29,  944.35,  797.34,  687.18,
              594.  ,  801.2 ,  821.79,  696.89,  572.3 ,  925.62,  788.07,
              934.75,  619.65,  805.29,  943.05,  882.55,  869.65,  991.82,
              724.01,  805.08,  838.01,  671.78, 1030.3 ,  785.95,  689.1 ,
              918.92,  618.3 ,  891.58,  899.  ,  776.67,  788.79,  891.47,
              828.73,  874.67,  796.43,  938.58,  588.55,  825.11,  828.27,
              733.4 ,  754.26,  836.85,  963.76,  751.67, 1021.7 ,  843.18,
              625.26,  970.49,  727.13,  762.51,  747.32,  608.57,  937.96,
              912.46,  805.94,  790.89,  743.31,  702.5 ,  901.4 ,  971.55,
              612.24,  780.07,  776.18,  804.51,  683.38,  994.11,  696.87,
              680.51,  830.42,  829.45,  721.85,  602.71,  845.69,  798.54,
              757.17,  800.49,  890.49,  850.44,  834.43,  573.53,  689.92,
              704.06,  692.2 ,  993.06,  948.83,  826.91,  700.09,  763.54,
              977.49,  772.38, 1034.7 ,  857.21,  773.6 ,  941.68,  780.37,
              782.98,  751.48,  832.19,  830.27,  889.23,  830.64,  746.92,
              804.48,  682.43,  558.6 ,  781.41,  959.49,  816.38,  860.01,
              702.16, 1115.3 ,  801.58,  867.54,  837.86,  792.1 , 1011.1 ,
              959.82,  829.68,  825.94,  699.84,  875.86,  835.  ,  889.76,
              746.42,  769.34,  686.17,  638.89,  968.98,  746.87,  951.02,
      ...
             1088.9 ,  953.25,  975.47,  739.3 ,  683.59,  762.3 ,  829.64,
              678.08, 1132.1 ,  927.19,  850.24,  920.17,  871.7 ,  923.98,
              724.22,  634.18,  778.44,  829.65,  653.52, 1084.7 ,  866.22,
             1104.4 ,  747.36,  802.2 ,  655.31,  745.85,  775.23,  777.17,
              974.98,  861.16,  542.66,  924.5 ,  583.76,  944.95,  708.19,
              711.34,  904.8 ,  832.63,  736.37,  755.84,  932.63,  877.49,
              835.05,  894.63,  738.95,  729.29,  728.25,  727.56,  856.33,
             1170.2 ,  677.33,  769.68,  917.17,  729.99,  906.88,  787.11,
              758.53,  891.55, 1056.2 ,  765.5 ,  935.82,  654.82, 1060.4 ,
              904.44,  754.59,  877.16,  705.33,  907.46, 1035.7 ,  778.12,
              952.24,  987.31,  710.32,  813.24,  915.55,  832.4 ,  690.63,
              924.76,  782.23, 1052.2 ,  915.47,  591.07,  580.45,  654.78,
              951.81,  849.29,  752.82,  657.91,  645.62,  818.47,  925.21,
              711.15,  977.86,  884.21,  814.93,  790.7 ,  832.93,  869.97,
              772.97,  861.66,  971.8 ,  774.22,  717.59, 1068.1 ,  902.09,
              779.6 ,  956.44,  890.78,  783.42,  863.74,  984.2 ,  689.48,
              896.96,  838.66,  788.  ,  923.13,  940.1 ,  993.93,  711.93,
              722.87,  805.08,  960.95,  762.37,  916.23,  758.47,  557.39,
              982.54,  898.36,  834.27,  822.31,  984.25,  813.69,  751.59,
             1047.6 ,  842.79,  799.76,  967.88])

    • ic_er
      (ic_er_dim_0)
      float64
      995.8 894.7 ... 1.123e+03 1.092e+03

      array([ 995.76,  894.69,  900.11,  831.04,  985.64,  862.56,  952.3 ,
             1012.3 , 1024.8 ,  874.49,  774.69, 1325.  ,  880.8 ,  718.72,
              736.77, 1092.4 ,  969.64,  853.82,  664.22, 1140.  ,  799.52,
             1009.5 ,  880.97,  806.59, 1130.5 , 1187.9 ,  942.49,  913.32,
              993.25,  941.41, 1144.4 ,  976.01, 1153.2 ,  963.54,  876.97,
              897.8 ,  819.38, 1070.1 , 1060.2 , 1054.  ,  833.05, 1148.3 ,
              938.83,  966.34, 1054.4 , 1098.1 ,  810.51, 1143.9 ,  929.43,
             1012.7 ,  994.12, 1026.5 , 1278.3 ,  906.09, 1162.7 ,  884.68,
              736.03, 1330.1 ,  878.61,  996.61, 1130.9 ,  810.97, 1100.5 ,
              950.35, 1304.  , 1115.9 ,  856.94,  752.47, 1197.4 , 1043.9 ,
              795.83,  946.68,  824.73,  831.37,  865.94, 1081.8 ,  916.09,
             1191.8 , 1059.3 , 1159.6 ,  921.97,  790.37,  866.67,  948.05,
              903.52, 1051.1 ,  923.99,  990.48,  845.59,  670.05,  964.66,
             1131.4 ,  848.99, 1075.7 , 1031.8 ,  993.83,  856.06,  980.45,
             1290.7 ,  862.07, 1185.4 , 1072.5 , 1023.6 ,  934.18,  902.04,
             1039.1 ,  909.56,  961.74,  963.76, 1107.1 ,  946.4 ,  917.13,
              822.85,  712.86,  714.01,  966.06, 1146.7 , 1120.  , 1048.7 ,
              984.69, 1225.5 , 1125.4 , 1004.8 , 1027.6 , 1015.1 , 1142.5 ,
              937.12,  879.88,  949.07,  875.49,  942.76,  957.79, 1028.2 ,
             1143.9 ,  979.5 ,  931.27,  722.43,  942.23,  851.91, 1082.1 ,
      ...
             1210.  , 1168.1 , 1210.7 ,  989.9 ,  935.29,  824.84, 1015.9 ,
             1121.  , 1174.5 , 1089.5 , 1116.4 , 1050.6 ,  910.56, 1069.6 ,
              993.27,  777.28, 1150.1 , 1038.7 ,  911.5 , 1275.4 ,  906.27,
             1416.7 ,  845.54, 1003.1 ,  891.1 ,  919.38,  854.6 ,  825.66,
             1122.1 , 1022.7 ,  835.91, 1084.8 ,  658.38, 1147.2 ,  757.77,
              794.42,  933.16,  965.96,  893.67,  912.42, 1123.3 , 1140.7 ,
              894.26, 1173.7 ,  762.06,  859.23,  803.77,  825.  , 1131.5 ,
             1417.8 ,  912.44,  766.79, 1034.7 ,  653.97, 1183.1 ,  865.21,
              704.74, 1012.5 , 1263.4 , 1057.5 , 1023.5 ,  622.99, 1304.7 ,
             1078.2 ,  972.28, 1050.2 ,  907.67,  907.25, 1017.8 ,  985.2 ,
              891.76, 1157.5 ,  910.15,  809.95, 1105.8 ,  993.83,  890.9 ,
             1175.6 ,  979.25, 1293.8 , 1026.3 ,  733.76,  777.04,  849.56,
             1138.7 , 1255.9 ,  845.15,  843.31,  931.4 , 1029.8 , 1165.4 ,
              871.97,  963.42, 1156.5 , 1150.3 , 1058.  , 1080.6 , 1080.4 ,
             1042.4 ,  973.06,  943.57,  865.61,  863.93, 1328.  , 1057.7 ,
              988.55, 1132.  ,  989.09,  685.23,  980.36, 1121.1 ,  840.17,
              966.41, 1175.4 , 1117.6 , 1265.8 , 1148.6 , 1067.5 ,  972.16,
             1169.7 ,  996.96,  965.23,  907.4 , 1067.7 , 1057.3 ,  650.59,
             1042.7 ,  988.99,  987.05, 1045.3 , 1246.5 , 1064.1 ,  869.78,
             1292.6 , 1041.3 , 1123.2 , 1091.7 ])

    • oc_ec
      (oc_ec_dim_0)
      float64
      664.0 624.3 526.9 ... 594.2 622.4

      array([664.008, 624.288, 526.872, 579.6  , 485.292, 477.864, 574.032,
             602.88 , 663.576, 638.664, 513.72 , 681.996, 488.052, 439.104,
             509.652, 531.792, 519.096, 435.588, 409.392, 674.028, 529.188,
             596.124, 500.784, 718.98 , 692.64 , 633.456, 584.388, 701.796,
             532.584, 648.384, 623.124, 506.784, 720.696, 642.576, 503.52 ,
             566.688, 537.24 , 591.828, 499.128, 589.32 , 632.616, 708.684,
             518.904, 521.532, 546.096, 584.244, 463.596, 565.524, 656.28 ,
             564.276, 501.384, 623.688, 671.916, 659.772, 647.748, 527.64 ,
             499.392, 712.848, 474.24 , 524.52 , 550.26 , 473.448, 666.216,
             585.408, 628.068, 578.172, 490.488, 473.796, 697.344, 586.788,
             458.232, 579.78 , 484.992, 523.872, 439.368, 648.528, 468.072,
             574.368, 648.804, 580.5  , 541.224, 518.736, 565.44 , 588.36 ,
             557.472, 589.644, 530.724, 496.536, 504.312, 493.68 , 599.184,
             589.488, 595.86 , 620.028, 728.496, 672.396, 522.744, 484.932,
             667.8  , 632.736, 662.784, 472.596, 579.276, 604.608, 546.744,
             583.536, 512.172, 470.628, 615.288, 682.608, 669.216, 528.156,
             562.332, 424.152, 505.056, 587.712, 664.548, 613.416, 636.456,
             654.984, 758.172, 578.304, 575.04 , 576.504, 560.592, 780.384,
             490.992, 586.86 , 582.9  , 567.3  , 607.98 , 499.668, 536.76 ,
             674.22 , 518.712, 598.872, 511.092, 561.876, 547.236, 557.58 ,
      ...
             729.12 , 691.02 , 647.256, 591.   , 585.54 , 486.444, 630.516,
             576.912, 688.512, 652.884, 651.06 , 603.876, 644.4  , 616.896,
             511.668, 427.152, 642.18 , 594.648, 456.348, 706.14 , 600.468,
             796.368, 486.396, 613.932, 439.536, 506.604, 510.324, 414.024,
             658.32 , 683.436, 380.7  , 517.08 , 423.552, 569.496, 441.888,
             603.576, 580.38 , 668.256, 531.972, 580.44 , 475.8  , 644.988,
             624.624, 664.44 , 399.588, 509.172, 434.292, 409.38 , 522.444,
             814.404, 496.116, 589.176, 683.58 , 500.172, 597.876, 531.948,
             449.4  , 660.204, 735.396, 658.956, 593.208, 386.94 , 661.584,
             610.68 , 394.104, 590.472, 597.36 , 566.676, 609.06 , 569.16 ,
             603.168, 614.148, 539.172, 477.228, 557.292, 513.18 , 479.64 ,
             709.668, 596.532, 743.748, 572.244, 466.632, 406.632, 488.16 ,
             632.676, 639.624, 494.016, 505.128, 514.92 , 539.1  , 647.304,
             471.324, 634.2  , 614.928, 585.468, 626.712, 562.776, 677.52 ,
             604.272, 658.152, 643.668, 474.516, 558.552, 773.664, 704.196,
             508.932, 720.432, 539.436, 588.96 , 642.936, 613.992, 522.372,
             603.24 , 596.772, 569.58 , 679.92 , 683.508, 602.784, 538.188,
             630.816, 524.808, 648.192, 568.248, 531.756, 609.78 , 446.16 ,
             598.536, 603.384, 481.284, 590.22 , 653.82 , 608.268, 569.376,
             846.432, 574.632, 594.24 , 622.416])

    • oc_er
      (oc_er_dim_0)
      float64
      606.7 583.8 485.7 ... 481.9 550.2

      array([606.72 , 583.788, 485.688, 510.264, 467.424, 445.248, 535.164,
             557.724, 542.112, 566.868, 456.936, 570.864, 405.636, 474.168,
             387.54 , 456.228, 535.656, 538.116, 348.888, 633.864, 343.548,
             497.748, 434.688, 467.52 , 612.792, 546.96 , 512.532, 567.72 ,
             528.252, 523.596, 596.028, 494.952, 660.264, 567.252, 438.192,
             505.908, 410.832, 523.5  , 568.38 , 519.096, 480.552, 696.876,
             479.748, 525.684, 491.664, 559.428, 355.692, 558.72 , 582.048,
             613.944, 517.752, 446.76 , 679.86 , 421.752, 584.82 , 451.968,
             439.788, 712.812, 479.844, 578.652, 409.14 , 489.132, 512.1  ,
             558.   , 555.012, 505.428, 471.924, 331.92 , 613.68 , 533.052,
             408.072, 489.732, 492.936, 560.292, 413.088, 494.928, 433.2  ,
             438.852, 420.168, 503.172, 471.684, 353.136, 515.94 , 569.664,
             509.472, 536.46 , 464.484, 477.768, 472.416, 460.572, 474.168,
             553.368, 497.496, 548.724, 586.32 , 574.2  , 526.404, 476.676,
             595.32 , 479.184, 634.272, 454.764, 546.144, 531.924, 479.784,
             536.76 , 517.584, 475.296, 565.344, 587.472, 521.268, 439.812,
             423.348, 284.004, 468.744, 376.212, 513.96 , 574.764, 574.32 ,
             484.908, 702.132, 530.832, 535.032, 507.372, 479.064, 645.84 ,
             483.972, 395.472, 499.956, 450.276, 516.864, 520.956, 516.468,
             556.884, 501.432, 513.684, 377.7  , 557.376, 499.932, 614.1  ,
      ...
             636.432, 585.42 , 673.332, 477.252, 486.312, 374.568, 479.928,
             534.396, 553.068, 605.112, 541.692, 519.756, 545.328, 544.38 ,
             526.176, 438.552, 593.088, 598.932, 394.584, 797.316, 503.556,
             685.86 , 449.736, 561.684, 422.76 , 477.42 , 528.276, 447.54 ,
             683.868, 577.956, 402.012, 464.784, 424.884, 516.432, 389.16 ,
             461.448, 444.012, 521.904, 427.008, 509.496, 449.28 , 589.416,
             458.88 , 532.44 , 506.16 , 449.652, 400.548, 357.744, 576.   ,
             730.716, 529.872, 478.98 , 603.684, 340.56 , 562.344, 542.436,
             422.58 , 530.268, 651.444, 427.308, 650.796, 440.52 , 707.772,
             561.96 , 406.44 , 497.232, 501.156, 480.876, 551.424, 563.46 ,
             565.5  , 619.176, 517.68 , 548.172, 547.332, 502.956, 437.364,
             598.656, 513.252, 638.868, 573.624, 369.432, 389.784, 398.76 ,
             568.212, 606.876, 491.568, 377.436, 509.028, 475.704, 656.064,
             400.752, 458.664, 673.512, 527.628, 446.244, 505.344, 549.672,
             655.86 , 561.024, 587.856, 487.74 , 412.308, 645.456, 624.312,
             436.968, 642.864, 572.004, 484.524, 481.344, 528.672, 440.82 ,
             495.936, 529.992, 496.212, 666.48 , 442.62 , 605.712, 531.216,
             414.684, 477.432, 506.34 , 556.272, 522.432, 556.464, 422.532,
             512.808, 524.892, 403.812, 546.756, 587.34 , 491.112, 504.78 ,
             661.932, 538.332, 481.872, 550.236])

  • Indexes: (4)
    • ic_ec_dim_0
      PandasIndex

      PandasIndex(Index([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,
             ...
             890, 891, 892, 893, 894, 895, 896, 897, 898, 899],
            dtype='int64', name='ic_ec_dim_0', length=900))

    • ic_er_dim_0
      PandasIndex

      PandasIndex(Index([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,
             ...
             890, 891, 892, 893, 894, 895, 896, 897, 898, 899],
            dtype='int64', name='ic_er_dim_0', length=900))

    • oc_ec_dim_0
      PandasIndex

      PandasIndex(Index([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,
             ...
             890, 891, 892, 893, 894, 895, 896, 897, 898, 899],
            dtype='int64', name='oc_ec_dim_0', length=900))

    • oc_er_dim_0
      PandasIndex

      PandasIndex(Index([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,
             ...
             890, 891, 892, 893, 894, 895, 896, 897, 898, 899],
            dtype='int64', name='oc_er_dim_0', length=900))

  • Attributes: (4)

    created_at :

    2024-06-23T20:34:06.404039+00:00

    arviz_version :

    0.18.0

    inference_library :

    pymc

    inference_library_version :

    5.15.1

  
  

In [27]:

idata_new_policy["predictions"]["p"].mean(dim=["chain", "draw", "obs"])

<xarray.DataArray 'p' (alts_probs: 5)> Size: 40B
array([0.05293203, 0.07065756, 0.66509458, 0.14841689, 0.06289894])
Coordinates:
  * alts_probs  (alts_probs) <U2 40B 'ec' 'er' 'gc' 'gr' 'hp'

xarray.DataArray

'p'

• alts_probs: 5

• 0.05293 0.07066 0.6651 0.1484 0.0629

  array([0.05293203, 0.07065756, 0.66509458, 0.14841689, 0.06289894])

• Coordinates: (1)
  • alts_probs
    (alts_probs)
    <U2
    'ec' 'er' 'gc' 'gr' 'hp'

    array(['ec', 'er', 'gc', 'gr', 'hp'], dtype='<U2')

• Indexes: (1)
  • alts_probs
    PandasIndex

    PandasIndex(Index(['ec', 'er', 'gc', 'gr', 'hp'], dtype='object', name='alts_probs'))

• Attributes: (0)

In [28]:

fig, ax = plt.subplots(1, figsize=(20, 10))
counts = wide_heating_df.groupby("depvar")["idcase"].count()
new_predictions = idata_new_policy["predictions"]["p"].mean(dim=["chain", "draw", "obs"]).values
ci_lb = idata_m3["posterior"]["p"].quantile(0.025, dim=["chain", "draw", "obs"])
ci_ub = idata_m3["posterior"]["p"].quantile(0.975, dim=["chain", "draw", "obs"])
ax.scatter(ci_lb, ["ec", "er", "gc", "gr", "hp"], color="k", s=2)
ax.scatter(ci_ub, ["ec", "er", "gc", "gr", "hp"], color="k", s=2)
ax.scatter(
    new_predictions,
    ["ec", "er", "gc", "gr", "hp"],
    color="green",
    label="New Policy Predicted Share",
    s=100,
)
ax.scatter(
    counts / counts.sum(),
    ["ec", "er", "gc", "gr", "hp"],
    label="Observed Shares",
    color="red",
    s=100,
)
ax.hlines(
    ["ec", "er", "gc", "gr", "hp"],
    ci_lb,
    ci_ub,
    label="Predicted 95% Credible Interval Old Policy",
    color="black",
)
ax.set_title("Predicted Market Shares under Old and New Pricing Policy", fontsize=20)
ax.set_xlabel("Market Share")
ax.legend();

Compare models

In [29]:

compare = az.compare({"m1": idata_m1, "m2": idata_m2, "m3": idata_m3, "m4": idata_m4})
compare

/Users/gzheng/projects/professional/likert_dualresponse/.venv/lib/python3.11/site-packages/arviz/stats/stats.py:789: UserWarning: Estimated shape parameter of Pareto distribution is greater than 0.7 for one or more samples. You should consider using a more robust model, this is because importance sampling is less likely to work well if the marginal posterior and LOO posterior are very different. This is more likely to happen with a non-robust model and highly influential observations.
  warnings.warn(

	rank	elpd_loo	p_loo	elpd_diff	weight	se	dse	warning	scale
m4	0	-1020.574185	9.992352	0.000000	7.557224e-01	28.247545	0.000000	False	log
m2	1	-1023.552608	5.038925	2.978423	2.305736e-01	27.799085	3.623620	False	log
m3	2	-1025.780290	9.345972	5.206104	1.972889e-13	28.163863	3.007272	True	log
m1	3	-1309.642978	1.228472	289.068792	1.370393e-02	12.918328	23.321276	False	log