Mixture of two beta geometric with beta shape > 1

Questions
1

I am trying to model a simple mixture of beta geometric distributions (following up from my previous question: Recovering parameters of long tailed Beta-Geometric data), but I get this error:

ValueError: Input dimension mis-match. (input[0].shape[0] = 2, input[1].shape[0] = 500)

my attempt this time:

import pymc3 as pm

model = pm.Model()

with model:    
    switch = pm.Uniform('switch', 0, 1)
    a1 = pm.Uniform('a1', 0.0001, 1)
    a2 = pm.Uniform('a2', 0.0001, 1)
    b1 = pm.Uniform('b1', 0.0001, 1)
    b2 = pm.Uniform('b2', 0.0001, 1)
    
    be1 = pm.Beta('be1',
                   alpha=a1,
                   beta=b1,
                   shape=population_size)

    be2 = pm.Beta('be2',
                   alpha=a2,
                   beta=b2,
                   shape=population_size)

    geo1 = pm.Geometric('geo1', p=(1-be1), shape=1)
    geo2 = pm.Geometric('geo2', p=(1-be2), shape=1)
    
    pm.Mixture('mixture',
       w=[switch, 1-switch], 
       comp_dists=[geo1.distribution, geo2.distribution],
       observed=geo
    )
    
    trace = pm.sample()
    
    trace = pm.sample()
    pm.traceplot(trace)

(the only way to make it work I found is to have 500 distributions in my mixture, but it makes little sense to me)

Thanks a lot in advance

1 Like
2

So I read this is related to an open bug. Not a big expert, but with some guidance is there anything I might contribute with to help solve this? @junpenglao

3

Hi @luca_giacomel, sorry about the delay response :sweat_smile: I tried your code some time ago and didnt manage to make it work, and it looks like a bug. I think in this particular example, it might be straightforward to write down the mixture logp directly. You can have a look of an example here:
http://docs.pymc.io/notebooks/gaussian-mixture-model-advi.html

4

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.

5

just tried, doesn’t change anything. I was calling “.dist” in the list of distributions so it’s pretty much equivalent

6

cool! :+1: I will try go under the hood and get my hands dirty with some logp stuff

7

Yeah I tried .dist but it doesnt work as well… :neutral_face: