Mixture of two beta geometric with beta shape > 1

[twiecki]

Shouldn’t you need to create the distributions with geo1 = pm.Geometric.dist(p=(1-be1), shape=1) See http://docs.pymc.io/api/distributions/mixture.html for an example.