Dirichlet with random shape

Questions
1

Hi,

I’m toying with a Normal MixtureModel. The weights are Dirichlet:

import pymc3 as pm
import numpy as np
import theano.tensor as tt
W = np.array([0.35, 0.4, 0.25])
with pm.Model() as model:
 w = pm.Dirichlet('w', np.ones_like(W)) 
 mu = pm.Normal('mu', 0., 10., shape=W.size)
 tau = pm.Gamma('tau', 1., 1., shape=W.size)
 x_obs = pm.NormalMixture('x_obs', w, mu, tau=tau, observed=x)

I am trying to change a bit the code so that instead of fixing the shape of the Dirichlet, the shape is a random uniform variable like this:

k = pm.DiscreteUniform(‘k’,lower=2,upper=5)
c = tt.ones_like(k)
w = pm.Dirichlet(‘w’, a=c)

but it is not working. The error is:

TypeError: Expected int elements in shape

I appreciate if you can help me how to go about this.
Many thanks

2

We do not allow dynamic shape in model, so currently this is impossible to do.
One potential work around is having a fix shape alpha with a shape is equal to the max shape (5 in this case), and dynamically assign zeros to some element. But again, inferencing these kind of model is also very tricky.

3

Thanks for the response.
Following your suggestion to dynamically assign zeros I have this attempt:

with pm.Model() as model:
 k = 5
 arr = []
 for j in range(k):
   arr.append(pm.Uniform('a'+str(j),lower=0,upper=0.99))
 w = pm.Dirichlet('w', a=tt.stack([arr[i] for i in range(k)]), shape=k)
 du = pm.DiscreteUniform('du',lower=0, upper=1,shape=k)
 mu = pm.Normal('mu', 0., 10., shape=k)
 tau = pm.Gamma('tau', 1., 1., shape=k)
 x_obs = pm.NormalMixture('x_obs', w, np.multiply(mu,du), tau=np.multiply(tau,du), observed=x)

Basically what I’m doing is to multiply by 0 or 1 the elements of w, mu, and tau. the 0's and 1's are those of the DiscreteUniform

It’s not working. The error is (only showing last fragment):

...........................................................................
/usr/local/lib/python2.7/dist-packages/pymc3/step_methods/hmc/base_hmc.py in astep(self=<pymc3.step_methods.hmc.nuts.NUTS object>, q0=array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
   0., 0.]))
144                 self.iter_count,
145                 None,
146                 None,
147             )
148             self._warnings.append(warning)
--> 149             raise SamplingError("Bad initial energy")
150 
151         adapt_step = self.tune and self.adapt_step_size
152         step_size = self.step_adapt.current(adapt_step)
153         self.step_size = step_size

SamplingError: Bad initial energy

I think the DiscreteUniform in particular lower=0 is generating the error.

If I change the bounds to something bigger than 0 it doesn’t throw the error but after sampling all du's are 1, defeating the purpose.

How can I implement the idea correctly?
Thanks

4

I have something that seemingly answers my previous question:

with pm.Model() as model:
	k = 5
	w = pm.Dirichlet('w', a=tt.stack([0.1,0.1,0.1,0.1,0.1]), shape=k)
	du = pm.DiscreteUniform('du',lower=0, upper=1,shape=k)
	idx_arr = []
	for i in range(k):
		if du[i]:
			idx_arr.append(i)

	mu = pm.Normal('mu', 0., 10., shape=k)
	tau = pm.Gamma('tau', 1., 1., shape=k)
	x_obs = pm.NormalMixture('x_obs', tt.stack([ w[h] for h in idx_arr ]), tt.stack([ mu[h] for h in idx_arr ]), tau=tt.stack([ tau[h] for h in idx_arr ]), observed=x)
	
with model:	
	trace = pm.sample(3000,n_init=10000, tune=500, random_seed=SEED)[1000:]

With respect to the first response, does this implements the suggestion of “dynamically assign zeros”?

5

Did this work for you @arraval - or did you find something else that worked? I’m facing the same issue…

6

@junpenglao Is there currently any support for dynamic shapes (i.e. using a tensorvariable as input to the shape argument)? I am trying to estimate the parameters of a Compound Poisson process with a Binomial distribution as the underlying distribution. The number of Binomial RVs follow a Poisson distribution, so I need a way to either dynamically set the shape argument or generate those RVs in a loop.

7

We are definitely better at supporting dynamic shape, but I am not sure if your use case is supported yet. Could you open a new post with the example you have in mind?

8

@junpenglao Sounds good. I have opened a new post here, thanks!