Problem with hierachical occupancy model

Questions
6

I would do it this way:
First reshape everything so within the pm.Model we are not dealing with dict:

index = []
W_r = []
y_r = []
for i in range(n):
    index.append(np.ones((W[i].shape[0], 1))*i)
    W_r.append(W[i])
    y_r.append(y[i][:, np.newaxis])
index = np.concatenate(index, axis=0).astype(int64)
W_r = np.concatenate(W_r, axis=0)
y_r = np.concatenate(y_r, axis=0)

Then I will use matrix multiplication instead of nested for loops:

with pm.Model() as m:
    beta = pm.Normal('beta', mu=0, sd=10, shape=(2, 1)) 
    alpha = pm.Normal('alpha', mu=0, sd=10, shape=(2, 1))
    x_b = pm.math.dot(X, beta)
    psi = pm.math.invlogit(x_b).flatten()
    
    w_a = pm.math.dot(W_r, alpha)
    d = pm.math.invlogit(w_a)
    
    # z = pm.Bernoulli('z', psi, shape=n)
    # prob = z[index] * d
    prob = psi[index] * d
    
    obs = pm.Bernoulli('y', p=prob, observed=y_r)
    trace = pm.sample(1000, tune=1000)

Note that I commented out the latent label part z = pm.Bernoulli('z', psi, shape=n), as the default value returns z being all zeros which gives error in the observed (it gives an error when you have p=0 but observed=1).

Then after inference, you can generate the latent label by sampling from a Bernoulli conditioned on the posterior:

with m:
    z = pm.Bernoulli('z', psi, shape=n)
    ppc = pm.sample_ppc(trace, vars=[z])

You can have a look at notebook here: https://github.com/junpenglao/Planet_Sakaar_Data_Science/blob/master/PyMC3QnA/discourse_1702.ipynb (which the second one use a mixture parameterization)

1 Like