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)
```

`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)
```