Showing that two models are equivalent

Questions version agnostic
1

I want to show that using a model with multiple observations and a Bernoulli distribution is equivalent to using the Binomial distribution. The real models are a bit more complex and involve mixtures, but I think it boils down to below example. The logp is probably different because of the constant in the binomial PMF so the two numbers shouldn’t be the same - but is there a way to use “something” that would show that they are equivalent?


with pm.Model() as m1:
    x = pm.Normal("x", mu=0, sigma=1)
    pm.Bernoulli("y", logit_p=x, observed=[0, 1, 0, 1, 1])

m1.logp().eval({'x': 0}) # returns -4.38467444

with pm.Model() as m2:
    x = pm.Normal("x", mu=0, sigma=1)
    pm.Binomial("y", logit_p=x, n=5, observed=3)

m2.logp().eval({'x': 0}) # returns -2.08208934
2

I don’t think there is a straightforward way to show that two models are equivalent.

I think you are right about the constant in the binomial PMF. When optimizing the calculation of logp, I am inclined to think that PyTensor decides to ignore any constant factors.

In your example, if you calculate the constant in the Binomial PMF (n over x), it is equal to (5 over 3) = 5! / (3! * 2!) = 10. And log(10) is approximately 2.3, which happens to be the difference between the values of logp that you get from the two models.

Maybe you can try to forse PyTensor to not do any optimizations to the computational graph and try to compare the logp again:

from pytensor.compile.mode import DebugMode

logp_1 = m1.logp()
logp_2 = m2.logp()

unoptimized_logp_1 = pytensor.function([], logp_1, mode=DebugMode())
unoptimized_logp_2 = pytensor.function([], logp_2, mode=DebugMode())

unoptimized_logp_1({'x': 0})
unoptimized_logp_2({'x': 0})

But, even with this, I am not sure it will include this constant or if it will work for every pair of equivalent models.

3

You can compile the logp and evaluate at 2 different points, the models should only differ by a constant.

Here is an example with a model that uses summary statistics and one that doesn’t

2 Likes
4

Thanks! Now that I think about it, another way to achieve the same is to look at gradients with dlogp() - kind of gets at the same thing.

5

In this case, you don’t need simulation—you can just do the derivation.

{N \choose \text{sum}(y)} \prod_{n=1}^N \text{bernoulli}(y_n \mid \theta) \\ = {N \choose \text{sum}(y)} \prod_{n=1}^N \theta^{y_n} \cdot (1 - \theta)^{1 - y_n} \\ = {N \choose \text{sum}(y)} \theta^{\text{sum}(y)} \cdot (1 - \theta)^{N - \text{sum}(y)} \\ = \text{binomial}(\text{sum}(y) \mid N, \theta).

So they only differ by a constant that does not depend on \theta.

I’m not sure that it’s enough to show the derivatives are the same—we need the likelihood to be the same for it to be guaranteed to produce the same result with HMC.

1 Like