How to access the logp function of a normal mixture?

Questions
1

Hi,

I’m having problem trying to access the logp function of a normal mixture.

This first attempt, let me think that I understand how to access the logp function of a basic random variable :

 with pm.Model() as model:
    test = pm.Normal('test',mu=0,sd=1)
    lp = model.logp({'test':12})

which gives lp =~ -73 (which sounds ok).

But I have a problem when trying to do the same stuff with a little more complex random variable, a normal mixture :

N = 5
with pm.Model() as model:
    w = pm.Dirichlet('w', a=np.ones(N))
    mu = pm.Normal('mu',mu=0,sd=1,shape=N)
    sd = pm.HalfNormal('sd',sd=1,shape=N)
    mix = pm.NormalMixture('mix',w=w,mu=mu, sd=sd) # building the mixture

    lp = model.logp({'mix':12}) # logp call

This gives me the following error : “TypeError: Missing required input: w_stickbreaking__” linked to the logp call…

I really do not understand why this doesn’t work: I see no functionnal difference between the two codes samples.

And the error message seem to talk about a missing feature w, which is explicitely given !! (Note : the code works without the logp call)

Any help will be appreciated.

2

What are you trying to do here? In general, you need to evaluate on a point that contains all the dependence (in this case would be all the free parameters in a model, which you can see a template by checking the model.test_point), but if you want to use the log_prob in a model, then you need to access the tensor (e.g., mix.logp(x))

3

I want to use the log_prob in my model (to use it later as a ‘potential’), so I gess that I have to access the tensor, something like :
mix.logp()

But I just can’t find the right syntax, all these attemps gives me errors :

lp = mix.logp(12)
lp = mix.logp([12])
lp =mix.logp(np.array[12])
lp = mix.logp({'mix':12})
lp = mix.logp({'mix':[12]})
lp = mix.logp({'mix':np.array([12])})

My real goal is to evaluate it at some given points created by a np.arange().

eval_points = np.arange(0,1,0.1)

and I’d like to have the logp of each individual point, not the sum of all logp.

How can I do that ?

4

If you are evaluating on some numpy array, why not just assign it to the observed?
Otherwise, you should do mix.distribution.logp(value)

5

Thanks very much, you suggestion :
mix.distribution.logp(value)
works, but I’m still curious about what does mix.logp then ? (it exists, so I guess it’s for some reasons…)

Which part of the documentation should I read to better understand this part ?

6

We does not expose these API, but basically mix.logp = mix.distribution.logp(mix) (recall that PyMC3 RandomVariables are theano tensors)

7

The logp is directly accessible from the RV name, and it seems to be different from the one accessible from the “distribution” sub class : see example below.

N = 5
with pm.Model() as model:
    w = pm.Dirichlet('w', a=np.ones(N))
    mu = pm.Uniform('mu',lower=0, upper=len(data),shape=N)
    sd = pm.HalfNormal('sd',sd=10,shape=N)
    mix = pm.NormalMixture('mix',w=w,mu=mu, sd=sd)

    i = np.arange(0,len(data))

    lp = mix.distribution.logp(i)   # this line works
    lp = mix.logp(i) # this line yield an error (see below)

The error message is :

TypeError: can’t turn [array([ 0, 1, 2, …, 1254, 1255, 1256])] and {} into a dict. cannot convert dictionary update sequence element #0 to a sequence

8

Correction: mix.logpt = mix.distribution.logp(mix) which is a tensor evaluated on itself
mix.logp is a compiled theano function that evaluated on all the dependency of mix (or all free parameters in model). Again you should not use it in model as it only evaluate once and will not add to the model logp (unless you wrap it in a pm.Potential)