Sample_ppc shape error with multi-index model

Questions
1

Hi,
I’m having a pesky problem with multi-indexing, which I never seem to be able to get right.
I have a multi-level beta-binomial model that models demand of certain product styles (10) over 52 weeks.
Each of the 10 styles has a different probability p for each of the 52 weeks, but I’d like to pool the styles across weeks which I try to do via the sd parameter of the priors on a and b of the beta on p.

WIth a lot of trial and error re: the indexes, I get the model to compile and efficiently sample - traceplots, gelman etc look good.
When I try to sample from the posterior I get the below error re: shapes. Any ideas? Any suggestions / examples on multi-indexing? I read through the radon and rugby examples already.

Sample Data:

(This only covers 5 styles over 13 weeks…hope that’s not confusing, I can repost the rest if that’d be more helpful)

style_demand = np.array([[1907, 1855, 1904, 1932, 1854, 1814, 1884, 1954, 1890, 1845, 1907,
        1913, 1856],
       [1528, 1531, 1572, 1593, 1567, 1542, 1591, 1586, 1578, 1573, 1605,
        1574, 1592],
       [1772, 1760, 1848, 1848, 1726, 1809, 1785, 1744, 1824, 1892, 1763,
        1748, 1804],
       [1065,  994, 1045, 1118, 1002, 1075, 1042, 1092, 1065, 1072, 1069,
        1046, 1073],
       [1605, 1650, 1705, 1753, 1694, 1664, 1663, 1719, 1697, 1623, 1751,
        1720, 1680]])

demand = np.array([ 9815,  9780,  9860, 10022,  9895,  9928,  9929, 10079,  9898,
       10030,  9980, 10034,  9870])

n_styles, n_weeks = style_demand.shape

Model:

with pm.Model() as model:
    
    # each of the 10 styles gets its own sds for and b, 
    # shared across all 52 weeks
    sd_a = pm.HalfNormal("sd_a", sd=1, shape=(n_styles, 1))
    sd_b = pm.HalfNormal("sd_b", sd=1, shape=(n_styles, 1))
    
    # each of the 10 styles has a different probability p 
    # in each of the 52 weeks
    a = pm.HalfNormal("a", sd=sd_a[style_idx,], shape=(n_styles, n_weeks))
    b = pm.HalfNormal("b", sd=sd_b[style_idx,], shape=(n_styles, n_weeks))
    
    style_probability = pm.Beta("style_probability", 
                                alpha=a, 
                                beta=b, 
                                shape=(n_styles, n_weeks))
    
    style_demand_obs = pm.Binomial("style_demand", 
                                   n=demand, 
                                   p=style_probability, 
                                   observed=style_demand 
                                  )

Shapes of the variables:

sd_a_log__ (10, 1)
sd_b_log__ (10, 1)
a_log__ (10, 52)
b_log__ (10, 52)
style_probability_logodds__ (10, 52)
style_demand (10, 52)

image
image.png892×920 38.7 KB

demand is the weekly overall sales number with shape (52,) and is a fixed array.

Error occurs here:

with model:
ppc = pm.sample_ppc(trace, vars=[style_demand_obs])

TypeError: Attempted to generate values with incompatible shapes:
        size: 1
        dist_shape: (10, 52)
        broadcast_shape: (10, 1)
2

Can confirm the error, but getting the shape to work is quite tricky. My suggestion is to flatten your predictors and observers into 10*52 vector - similar hierarchical model are much easier to handled and there are quite some example around.

3

FYI, the error is the broadcasting in HalfNormal:

a.random(point=trace[0])

will generate the same error

4

Oh I see, and then just index into the relevant rows for each style. That makes sense, I’ll give that a shot. I think I often think about the hierarchical set up too literally.
Thanks as always!

5

No problem! I filed a bug on Github: https://github.com/pymc-devs/pymc3/issues/3147

I like the 1d vector set up better, because it is also easier if you have missing data.

6

This appears to work beautifully. Thanks again!

with pm.Model() as model:
    
    # each of the styles gets its own sds for and b, 
    # shared across all weeks
    sd_a = pm.HalfNormal("sd_a", sd=2, shape=n_styles)
    sd_b = pm.HalfNormal("sd_b", sd=2, shape=n_styles)
    
    # each of the styles has a different probability p 
    # in each of the weeks
    a = pm.HalfNormal("a", sd=sd_a[style_idx], shape=(len(style_idx),))
    b = pm.HalfNormal("b", sd=sd_b[style_idx], shape=(len(style_idx),))
    
    style_probability = pm.Beta("style_probability", 
                                alpha=a[style_idx], 
                                beta=b[style_idx], shape=(len(style_idx),))
    
    style_demand_obs = pm.Binomial("style_demand", 
                                   n=demand, 
                                   p=style_probability, 
                                   observed=style_demand 
                                  )

Where all inputs are now flattened to 1D vectors.

1 Like