Beta-Binomial conjugate prior -- pm.Binomial buggy results...?

I am not sure where did I make a typo. A beta-binomial distribution or model refers to this kind of model where each observation is modelled with a binomial likelihood whose p parameter is in turn given a unique beta prior that has shared alpha/beta hyper priors.

This:

with pm.Model() as m1:
  alpha = pm.HalfNormal('alpha',sigma=1)
  beta = pm.HalfNormal('beta',sigma=1)
  theta = pm.Beta('theta',alpha, beta, shape=len(W))
  lik = pm.Binomial('likelihood',n=N,p=theta,observed=W)

And this are therefore equivalent:

with pm.Model() as m1:
  alpha = pm.HalfNormal('alpha',sigma=1)
  beta = pm.HalfNormal('beta',sigma=1)
  lik = pm.BetaBinomial('likelihood',n=N,alpha=alpha, beta=beta,observed=W)

Except theta is marginalized out in the later

2 Likes