Adding constraint to the parameter

Questions
1

Hi,

I am trying to make a scenario where three uniformly distributed parameters need to be constrained.
p1 = pm.Uniform(‘p1’,lower = 0, upper = 1)
p2 = pm.Uniform(‘p2’,lower = 0, upper = 1)
p3 = pm.Uniform(‘p3’,lower = 0, upper = 1)

I want the constraint to be p1+p2+p3 = 1. How can I do it?

p = np.array([p1,p2,p3])
p_potential = pm.Potential(‘p_min_potential’, tt.switch(tt.sum§ = 1., 0, -np.inf))

I tried using this but it gives error because in my knowledge the data type of p is something which is not tensor type.

Help much appreciated.
Thanks

2

§ is (p) in my question.

3

What you are describing is a simplex, which is usually modelled using a Dirichlet distribution:

p = pm.Dirichlet('p', a=np.ones(3))
4

Thanks @junpenglao
I was wondering if there is a way without using dirichlet and adding a constraint to my problem using pm.Potential(‘p_min_potential’, tt.switch(tt.sum(p) = 1., 0, -np.inf)) .

Can we even do this kind of operation on parameters with some distribution ?

5

You can use potential sure but sampling will be really inefficient.

pm.Potential(‘p_min_potential’, tt.switch(tt.eq(tt.sum(p), 1.), 0, -np.inf))

You can also normalized it p = p/p.sum() but Dirichlet is the best solution.

6

@junpenglao what do you mean by p here ?
Because if I use
p = np.array([p1,p2,p3]) I am going to get same error.

7

Try p = tt.stack([p1, p2, p3])

8

@junpenglao Thanks for help. I got the code running now but results are not making any sense. I am getting the distribution of p with its mean sum much greater than 1.

9

The problem with p = p/p.sum() is that the solution will not be unique: [.1, .1, .1] == [.2, .2, .2] under that setup.
Setting one of the p’s to 1. would make it unique :

p_ = tt.concatenate([pm.HalfFlat('p_', shape=2), np.array([1.])])
p = pm.Deterministic('p', p_/p_.sum())

Sorry, not the most beautiful code.
You could also use a masked array with the observed parameter of pm.HalfFlat.

1 Like
10

Your solution worked but I did not understand the part about why p = p/p.sum() is not unique. Could you please elaborate on that?

Thanks

11

With b = a/a.sum(), you can multiply any solution of the vector a by any scalar without affecting b. ie, b = [1,1,1]/(1+1+1) yields the same end value as b = [.5,.5,.5]/(.5+.5+.5). This means that the inference will drift endlessly between all and any homothetic solution for a.

Hope that helps?

12

I understood what you are trying to tell here but how does making one of the p’s 1 is making any difference. I can’t see that.

13

Yep - which is why you should use Dirichlet, because for a n simplex it is defined on n-1 parameters that use a similar treatment.

14

In a n element simplex, since the sum is define as 1, when you specified the first n-1 element, the last element is already know as you can do p_n = 1-p_1-p_2-...-p_{n-1}, which means you only need to specify the first n-1 element otherwise the model is over-parameterised.

15

Any value would have worked BTW, it’s just that 1 is neat.

Say that the correct solution of p is [.3,.2,.5], if you fix the last value to 1, the other values can only be [.6, .4]

16

So what you are suggesting is that even if I have more than 3 p’s , to have same constraint it is a good option to use Dirchlet rather than uniform for the prior?

17

Absolutely. Dirichlet is much better tested with a stable internal transformation to make sure model is identifiable (ie no need to worry about specifying n-1 free variables)

18

@junpenglao @gBokiau Thanks a lot for your guidance.
I really appreciate it.

1 Like