Folks I am facing similar issue, my code (@cluhmann and @ckrapu

```
import pymc3 as pm
import arviz as az
# ri~ Binomial(pi, ni)
#
# logit(pi) = α0 + α1x1i + α2x2i + α12x1ix2i + bi
#
# bi~ Normal(0, τ)
import numpy as np, pymc3 as pm, theano.tensor as T
# germinated seeds
r = np.array([10, 23, 23, 26, 17, 5, 53, 55, 32, 46, 10, 8, 10, 8, 23, 0, 3, 22, 15, 32, 3])
# total seeds
n = np.array([39, 62, 81, 51, 39, 6, 74, 72, 51, 79, 13, 16, 30, 28, 45, 4, 12, 41, 30, 51, 7])
# seed type
x1 = np.array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1])
# root type
x2 = np.array([0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1])
# number of plates
N = x1.shape[0]
import math
with pm.Model() as m:
### hyperpriors
sigma = pm.Uniform('sigma',lower=0.,upper=10.)
tau =1/pow(sigma,2)
### parameters
# fixed effects
alpha_0 = pm.Normal('alpha_0', mu=0.,tau=1e-6)
alpha_1 = pm.Normal('alpha_1', mu=0.,tau=1e-6)
alpha_2 = pm.Normal('alpha_2', mu=0.,tau=1e-6)
alpha_12 = pm.Normal('alpha_12', mu=0.,tau=1e-6)
# random effect
b = pm.Normal('b', 0., tau, shape=(N,))
# expected parameter
logit_p = (alpha_0 + alpha_1 * x1 + alpha_2 * x2 + alpha_12 * x1 * x2 + b)
p = T.exp(logit_p) / (1 + T.exp(logit_p))
### likelihood
obs = pm.Binomial('obs', n, p, observed=r)
n = 10000
with m:
start = pm.find_MAP({'sigma':8., 'alpha_0': 0., 'alpha_1': 0., 'alpha_2': 0., 'alpha_12': 0.})
step = pm.HamiltonianMC(scaling=start)
ptrace = pm.sample(n, step, start, progressbar=False)
start = pm.find_MAP({'sigma':1., 'alpha_0': 0., 'alpha_1': 0., 'alpha_2': 0., 'alpha_12': 0.,'b': np.zeros(N)})
step = pm.HamiltonianMC(scaling=start)
ptrace = pm.sample(n, step, start, progressbar=False)
burn = 1000
pm.summary(ptrace[burn:])
```

It is nowhere close to actual params. Any help appreciated