Basically I’m trying to run a Beta-Binomial regression with categorical input features that I factorize & add to coords, and some purely numerical features that are multiplied with beta.

I’ve been fiddling with this model for the past couple of weeks & every time I run it the tuned samples don’t diverge but when the first 500 iterations are done it starts diverging for the majority of the non-tuned samples.

I’ve tried playing around with the NUTS sampler by increasing the target_accept to no avail. I would like to re-parameterize if that will help but I don’t know which part of my code might be giving the issue.

You will find the piece of code below (apologies this is a class method on an object I haven’t included). Can anyone give me a hint where I should look? As input I put the following:

- self.draws = 500
- 7x numerical features that are multiplied with beta after being standardized
- 16x categorical features that are turned into "_alpha"s. The largest cardinality of which is 39 for one feature

```
# Factorize categorical features
cat_idxs = dict()
for feat in self.cat_feats:
cat_idxs[feat], cat_vals = pd.factorize(X[feat])
if len(self.coords) == 0:
self.coords = {"$obs_id$": np.arange(len(cat_idxs[feat])), "$num_feats$": self.num_feats, feat: cat_vals}
else:
self.coords[feat] = cat_vals
# Construct model
with pm.Model(coords=self.coords) as model:
# Construct Numerical Beta Priors
beta, x_ = None, None
if self.num_feats:
beta = pm.Normal(name="beta", mu=0, sigma=20, dims="$num_feats$")
# Construct Data Containers
x_ = pm.Data("X", X[self.num_feats])
# Construct Categorical Alpha Priors
cat_alphas = dict()
for feat, cat_vals in self.coords.items():
if feat not in ["$obs_id$", "$num_feats$"]:
cat_idxs[feat] = pm.Data(f"{feat}_idx", cat_idxs[feat], dims="$obs_id$")
cat_alphas[feat] = pm.Normal(name=f"{feat}_alpha", mu=0, sigma=20, dims=feat)
# Construct "variance" prior(s)
# When kappa is large, the beta is less variable
kappa = pm.HalfNormal("kappa", sigma=20)
# Construct Deterministic Target Mean
if not self.cat_feats:
alpha = pm.Normal(name="alpha", mu=0, sigma=20)
mu = pymc3_math.dot(x_, beta) + alpha
else:
if self.num_feats:
mu = pymc3_math.dot(x_, beta)
else:
mu = 0
for feat in self.cat_feats:
# Add the categorical alphas to mu
mu += cat_alphas[feat][cat_idxs[feat]]
# Define mean with sigmoid transform (-inf, inf) -> (0, 1)
mu = pm.Deterministic(name="mu", var=pm.math.invlogit(mu))
# Pass this sigmoid mean to beta distribution
ps = pm.Beta('ps', mu*kappa, (1-mu)*kappa, dims="$obs_id$")
# Pass beta to binomial likelihood w observations
_ = pm.Binomial('target', n=10, p=ps, observed=y)
logging.info("Start Sampling...")
step = pm.NUTS(
target_accept=.95
)
self.trace = pm.sample(draws=self.draws, step=step, tune=self.draws, n_init=200_000)
self.burned_trace = self.trace[int(0.2*self.draws):]
```