Hi all,

I am new to PyMC3. I am having trouble implementing the occupancy model shown in the images above. The alpha and beta each have a prior distribution that follows a multivariate normal distribution.

```
with pm.Model() as _sampler:
beta = pm.Normal('beta', mu=0, sd=1000, shape=len(HYPERS['b_mu']))
alpha = pm.Normal('alpha', mu=0, sd=1000, shape=len(HYPERS['a_mu']))
z = []
y_ = []
for i in range(X.shape[0]):
x_b = pm.math.dot(X[i], beta)
psi = pm.math.invlogit(x_b)
z.append(pm.Bernoulli('z'+str(i), psi))
ydata = y[i]
yy = []
for j in range(W[i].shape[0]):
w_a = pm.math.dot(W[i][j], alpha)
d = pm.math.invlogit(w_a)
prob = z[i] * d
yy.append(pm.Bernoulli('yy'+str(i)+','+str(j), p=prob, observed=ydata))
y_.append(yy)
z = tt.stack(z)
y_ = tt.stack(y_)
with _sampler:
# obtain starting values via MAP
start = pm.find_MAP(model=_sampler)
# instantiate sampler
step = pm.Metropolis()
_trace = pm.sample(5000, step=step, start=start)
pm.traceplot(_trace[2500:])
```

The X is a N by p matrix that has the covariate information regarding the beta parameter and W is a N by K by q 3 dimensional array that contains the covariate information regarding the alpha parameter. y is stored in python as a dictionary in Python with its keys being the index (1 to N). The W array is also stored as a dictionary where each element of the keys is a 2d numpy array. Also in the code I used a different link function than in the pictures (logit instead of probit).

The problem is when I run this code this is the output I get from the console:

logp = -inf, ||grad|| = 12.207: 100%|##############################################| 3/3 [00:00<00:00, 28.29it/s]

{‘beta’: array([ -8.19231921e-01, -8.18577127e-17]), ‘alpha’: array([ 0.57346234, 0. ]), ‘z0’: array(0), ‘z1’: array(0), ‘z2’: array(0), ‘z3’: array(0), ‘z4’: array(0), ‘z5’: array(0), ‘z6’: array(0), ‘z7’: array(0), ‘z8’: array(0), ‘z9’: array(0), ‘z10’: array(0), ‘z11’: array(0), ‘z12’: array(0), ‘z13’: array(0), ‘z14’: array(0), ‘z15’: array(0), ‘z16’: array(0), ‘z17’: array(0), ‘z18’: array(0), ‘z19’: array(0)}

After that nothing happens and it hangs forever. I can implement this algorithm by directly working out the posterior full conditionals and build a class where I can iteratively sample using Gibbs, but I want to make it work using PyMC3 but it seems like I am doing something wrong. Please help!