pm.Categorical with two p's, not giving same results as pm.Bernoulli

Questions
1

I could be doing something stupid, or I could just not be understanding the math correctly between the two … but my understanding is that if you give the Categorical distribution two probabilities (p(x=0) = 1-p_success, p(x=1) = p_success), that should be exactly the same thing as the Bernoulli distribution. So why aren’t traces from two otherwise identical models (and same random_seed) not giving equal results??

obvs_ind5 = [1,1,1,1,1,1,1,0,0,0] # p_true_ind5 = 0.7
ind = 5

with pm.Model() as model1:
    z = pm.Normal('z',mu=0,sd=1,shape=(10,))
    p1 = pm.math.invlogit(z)
    for i in range(len(obvs_ind5)):
        pm.Bernoulli('y_%d'%(i),p=p1[ind],observed=obvs_ind5[i])
    trace = pm.sample(progressbar=False,random_seed=123)
Z1 = trace['z']

with pm.Model() as model2:
    z = pm.Normal('z',mu=0,sd=1,shape=(10,))
    p2 = pm.math.invlogit(z)
    for i in range(len(obvs_ind5)):
        pm.Bernoulli('y_%d'%(i),p=p2[ind],observed=obvs_ind5[i]) # TRUE
#         pm.Categorical('y_%d'%(i),p=[1-p2[ind],p2[ind]], observed=obvs_ind5[i]) # FALSE
#         pm.Categorical('y_%d'%(i),p=[1-p2[ind],p2[ind]], observed=1-obvs_ind5[i]) # FALSE
#         pm.Categorical('y_%d'%(i),p=[p2[ind],1-p2[ind]], observed=obvs_ind5[i]) # FALSE
#         pm.Categorical('y_%d'%(i),p=[p2[ind],1-p2[ind]], observed=1-obvs_ind5[i]) # FALSE
    trace = pm.sample(progressbar=False,random_seed=123)
Z2 = trace['z']

print((Z1 == Z2).all())

I tried all combinations I could think of for the parameters and observation, but it never works. Obviously setting both to Bernoulli gives the True result. Any help?

2

The two model is the same, I simplify your model a bit:

with pm.Model() as model1:
    z = pm.Normal('z', mu=0, sd=1)
    p_ind5 = pm.Deterministic('p', pm.math.invlogit(z))
    obs = pm.Bernoulli('y', p=p_ind5, observed=obvs_ind5)

with pm.Model() as model2:
    z = pm.Normal('z', mu=0, sd=1)
    p_ind5 = pm.Deterministic('p', pm.math.invlogit(z))
    obs = pm.Categorical('y', p=[1-p_ind5, p_ind5], observed=obvs_ind5)

Sample from it would give you the same posterior distribution of z. Of course, you cannot expect to get identical trace, a better way to check the model is identical is to check the logp directly:

for i in np.linspace(-1, 1, 10):
    np.testing.assert_almost_equal(model1.logp({"z": i}), model2.logp({"z": i}))
3

Hi @junpenglao, thanks for reply! Two question for you:

1 ) The main thing you did differently in model is to use pm.Deterministic on the invlogit. Is that necessary and/or recommended? I didn’t use it in the larger version of code I’ve built before (using pm.Bernoulli) and the code worked. I just gave the pm.Bernoulli the “p” directly, without any pm.Deterministic “wrapper”. Do you recommend it?

2 ) Your statement “you cannot expect to get identical trace” … but why not?? IF the random_seed of the pm.trace is the same, and it truly is the same model, then why should it not produce the same trace?? It DOES do this when model1 and model2 using exactly the same pm.Bernoulli for the obvs.

Thank you again!!

4

It’s not necessary - was just to have it in the trace.

The two logp function is written differently, so when theano compile the model logp it might optimized the graph differently.

1 Like
5

Thank you!! @junpenglao