Add bounds to gaussian unmixture model

Questions v5
1

I have a dataset that consists of two truncated normal distributions that I’m trying to unmix. I know the mu/sigma of one of the two and I’m trying to extract the weights, as well as the mu/sigma of the unknown gaussian. The follow code works most of the time:

with pm.Model() as model:
                w = pm.Dirichlet('w', a=np.array([1, 1]))  # 2 mixture weights
                mu1 = pm.Normal("mu1", mu=np.mean(x), sigma=np.std(x))
                sig1 = pm.HalfNormal('sig1', sigma=np.std(x))
                
                tn1 = pm.TruncatedNormal.dist(mu=mu1,
                                            sigma=sig1,
                                            lower=low_cutoff,
                                            upper=high_cutoff)
                
                tn2 = pm.TruncatedNormal.dist(mu=known_mean,
                                            sigma=known_sigma,
                                            lower=low_cutoff,
                                            upper=high_cutoff)
                
                like = pm.Mixture(name="like",
                                w=w,
                                comp_dists=[tn1, tn2],
                                observed = x)

However sometimes it converges on some obviously wrong results. I know things about the distribution such as the bound on the weights (no distribution should be less than 5% of the mixture, sig1 should be bounded between 0.5 and 2), but I can’t quite figure out how to add bounds like this.

I tried:

pm.Bound(pm.HalfNormal, lower=0.5, upper=2.0)('x', mu=np.mean(x), sigma=np.std(x))

and I get an error about dist:

TypeError: __new__() missing 1 required positional argument: 'dist'

It seems like this package would have a nice way to give the model this range before hand, but I just can’t quite figure it out in the documentation/googling around.

2

I suggest you use pm.Truncated instead of the old pm.Bound: Truncated — PyMC dev documentation.

There is an examples in the docs of the API. Bound is actually similar if you really insist on it

3

Cor your weights you can give a slightly larger alpha prior like [10, 10], that is a bit more compatible with your knowledge that no distribution can be less than 5. Instead it’s a bit more concentrated towards the 50-50 split

Otherwise just scale the weights, like w = 0.05 + w*0.9 or something more gradual that goes in the same direction

4

So I tried this about an hour ago and maybe I just didn’t quite get it right? But it gave me a solution outside of the bounds I have given:

sig1 = pm.Truncated('sig1', pm.Normal.dist(mu=np.mean(x), sigma=np.std(x)), lower=0.50, upper=2.0)

Gave me values for sig1 that were all at 3.5ish (and were not good fits with the underlying data).

Did I mess that up somehow?

5

It’s not possible to give you values outside of the range