Why can't I use a Bernoulli as a likelihood variable in a hierarchical model in PyMC3?


#1

#2

The shape error you are having is because you broadcasted the alpha and beta of the pm.Beta too early. In

    # Parameters for coins
    theta = pm.Beta('theta',
                     omega_c[mint_idx]*(kappa_c[mint_idx]-2)+1,
                    (1-omega_c[mint_idx])*(kappa_c[mint_idx]-2)+1,
                     shape = num_coins)

You want the theta has the shape of num_coins, and then broadcasted it during observed. Because the latent variables (theta of each coin) is elementary/unique.
So in you first example, you should do something like:

grouped = observations.groupby(['mint', 'coin']).agg({'outcome': [np.sum, np.size]}).reset_index()
grouped.columns = ['mint', 'coin', 'heads', 'total']

num_mints = grouped['mint'].nunique()
mint_idx = grouped['mint']
num_coins = observations['coin'].nunique()
coin_idx = observations['coin']

with pm.Model() as hierarchical_model2:
    # Hyper parameters
    omega = pm.Beta('omega', 1, 1)
    kappa_minus2 = pm.Gamma('kappa_minus2', 0.01, 0.01)
    kappa = pm.Deterministic('kappa', kappa_minus2 + 2)

    # Parameters for mints
    omega_c = pm.Beta('omega_c',
                       omega*(kappa-2)+1, (1-omega)*(kappa-2)+1,
                       shape = num_mints)    
    kappa_c_minus2 = pm.Gamma('kappa_c_minus2',
                              0.01, 0.01,
                              shape = num_mints)
    kappa_c = pm.Deterministic('kappa_c', kappa_c_minus2 + 2)

    # Parameters for coins
    theta = pm.Beta('theta',
                     omega_c[mint_idx]*(kappa_c[mint_idx]-2)+1,
                    (1-omega_c[mint_idx])*(kappa_c[mint_idx]-2)+1,
                     shape = num_coins)

    y2 = pm.Bernoulli('y2', p=theta[coin_idx], observed=observations['outcome'])

#4

Ah, it was the shape of mint_idx. Thank you so much!