Hi,

New to PYMC here.

I’ve been trying to follow A Primer on Bayesian Methods for Multilevel Modeling — PyMC3 3.11.5 documentation

I see that to predict on test data, I need to do:

```
prediction_coords = {"obs_id": ["ST LOUIS", "KANABEC"]}
with contextual_effect:
pm.set_data({"county_idx": np.array([69, 31]), "floor_idx": np.array([1, 1])})
stl_pred = pm.fast_sample_posterior_predictive(
contextual_effect_idata.posterior, random_seed=RANDOM_SEED
)
az.from_pymc3_predictions(
stl_pred, idata_orig=contextual_effect_idata, inplace=True, coords=prediction_coords
)
```

However, how do I do a prediction if my county_idx is new and never seen before in the training data?

I tried to a pass a new county_idx and get the following error:

```
IndexError: index 10 is out of bounds for axis 0 with size 10
Apply node that caused the error: AdvancedSubtensor1(intercept_county, county_id_index)
Toposort index: 0
Inputs types: [TensorType(float64, (None,)), TensorType(int32, (None,))]
Inputs shapes: [(10,), (10,)]
Inputs strides: [(8,), (4,)]
Inputs values: ['not shown', 'not shown']
Outputs clients: [[normal_rv{0, (0, 0), floatX, True}(RandomGeneratorSharedVariable(<Generator(PCG64) at 0x15C2F2F20>), TensorConstant{[]}, TensorConstant{11}, AdvancedSubtensor1.0, sigma)]]
HINT: Use the Aesara flag `exception_verbosity=high` for a debug print-out and storage map footprint of this Apply node.
```

I saw How do we predict on new unseen groups in a hierarchical model in PyMC3? - #2 by lucianopaz but not sure how to apply that to this simple model:

```
with pm.Model(coords=coords) as varying_intercept:
floor_idx = pm.Data("floor_idx", floor, dims="obs_id")
county_idx = pm.Data("county_idx", county, dims="obs_id")
# Hyperpriors:
a = pm.Normal("a", mu=0.0, sigma=10.0)
sigma_a = pm.Exponential("sigma_a", 1.0)
# Varying intercepts:
a_county = pm.Normal("a_county", mu=a, sigma=sigma_a, dims="County")
# Common slope:
b = pm.Normal("b", mu=0.0, sigma=10.0)
# Expected value per county:
theta = a_county[county_idx] + b * floor_idx
# Model error:
sigma = pm.Exponential("sigma", 1.0)
y = pm.Normal("y", theta, sigma=sigma, observed=log_radon, dims="obs_id")
pm.model_to_graphviz(varying_intercept)
```

One idea I had is to train with an additional fake index (which won’t show up in observations) and assign any new counties to that index. That way it should still get the benefit of the partial pooling and have a prediction (albeit highly uncertain).

Any help is appreciated!