# Fitting Gaussian mixture model with constraints (eg. mu1<mu2) in Python

My question is similar to this one, but while the OP there has constrains such as mu1 being <=0 and mu2 being >=0, my constraints are following:

• It's a three component mixture model.
• mu1 < mu2 < m3
• mu3>=2 x mu1
• w1=w3 and w2<w1 and w2<w3

The solution that is posted is one using autograd python module. We define negative_log_likelihood as the function we want to minimize:

def negative_log_likelihood(params,data):
w1,m1,s1,w2,m2,s2,w3,m3,s3=params
return(-np.log(w1*norm.pdf(data,m1,s1)+w2*norm.pdf(data,m2,s3)+w3*norm.pdf(data,m3,s3)).sum())


We import necessary modules (minimize and value_and_grad):

from autograd import value_and_grad
from scipy.optimize import minimize


and then comes the actual function where I'm stumped:

results = minimize(
x0 = np.array([1, 0, -1, 0, 0.5]), # initial value
args=(obs,),
jac=True,
bounds=(
(0, None),    # mu1 (you mentioned the constraints on the means)
(None, None), # log_sigma1 is unbounded
(None, 0),    # mu2 (you mentioned the constraints on the means)
(None, None), # log_sigma2 is unbounded
(0, 1)        # the weight param should be between 0 and 1
)


)

This part of the code I C/P from the original question, the OP has 5 parameters (well 6) and I have 9 (3 for each component). So autograd value_and_grad works on negative likelihood function, then scipy.optimize module minimize, minimizes it. X0 is our initial guess, and args is our data array. I am stumped on how to write the bounds section, and if that will even work.

EDIT#1:

if I have m1<m22m1 Reformating: (1) 0<m2-m1 (2) 0<m3-2m1

(3) 0<m3-m2 (4) 0<m3-2m1

Combining (1) and (2), ditto for 3 and 4:

0<m3-3m1+m2 REFORMATING : -3m1+m2+m3 0<2m3-m2-2m1 REFORMATING : -2m1-m2+2m3

and I get [-3,1,1] as first row, [-2,-1,2]

EDIT #2

Code I'm making is basically code by @Cam.Davidson.Pilon with slight modifications (modules imported, not shown):

def mixgdata_make(pi1=0.45,pi2=0.1,pi3=0.45,m1=200,m2=300,m3=400,s1=50,s2=100,s3=50,size=1000):
#this works
import numpy as np
result=np.concatenate((np.random.normal(loc=m1,scale=s1,size=int(size*pi1)),np.random.normal(loc=m2,scale=s2,size=int(size*pi2)),np.random.normal(loc=m3,scale=s3,size=int(size*pi3))),axis=0)
np.random.shuffle(result)
return(result)


This gives me NumPy Object Array with type float 64, with size (1000,). So far so good. I call the data with:

data=mixgdata_make()


I define my negative log likelihood function similar to OP:

def negative_log_likelihood(params,data):
m1,m2,m3,w1,w2,w3,s1,s2,s3=params
return(-np.log(w1*norm.pdf(data,m1,s1)+w2*norm.pdf(data,m2,s3)+w3*norm.pdf(data,m3,s3)).sum())


This also works, as I've tested it with bad parameters (way off) and with good parameters (pretty close).

But here is the problem:

results=minimize(fun=value_and_grad(negative_log_likelihood),
x0=np.array([1,0,-1, 0,0.5,1, 2,4,5]),
args=(data,),
)


I get ValueError. But if I set jac=True, it works! In a sense that the code runs. But optimization itself has status as 2, success as False, x as basically what is inputed and message of desired error not necessarily achieved due to precision loss.

So while I implemented it (without Constraints, but one step at a time), it doesn't actually work.

Maybe the error is in supplied x0? My first idea will be to try to do Gaussian Mixture Model in scipy on original data to get "estimates" of all nine parameters, then feed those parameters in x0. I'll edit this post with an answer.

• Can you combine inequalities like this? I'm not so sure... – Cam.Davidson.Pilon Feb 22 at 18:39
• @Cam.Davidson.Pilon I don't see why not, but I don't really know enough to say for sure. Hence coming here. Also, if you don't mind, I've taken the liberty to take code you originally wrote for another person and see if it will work for my case. However, it's giving me the error ValueError: setting an array element with a sequence.. I'll edit my original post to show the code. Could I trouble you for additional help? – ivan199415 Feb 23 at 12:03
• 0<m2-m1 & 0<m3-2m1 implies 0<m3-3m1+m2 but the converse isn't true. (Simpler example: 0 < x & 0 < y implies 0 < x+y but 0 < x+y doesn't mean that both 0 < x and 0 < y. – Cam.Davidson.Pilon Feb 23 at 23:41
• But if I set jac=True, it works! right, I'm not surprised. By saying jac=True, you are saying to minimize that the function given returns two outputs (look up the docs on jac=True). – Cam.Davidson.Pilon Feb 23 at 23:43
• Convergence problems (i.e status=2) are tricky to debug, and solutions can be bespoke to your model. Try non-default optimization algorithms that minimize provides (see method kwarg) – Cam.Davidson.Pilon Feb 23 at 23:44

Looking at scipy's minimize docs, we can use LinearConstraint in the constraints kwarg (instead of bounds).

First, we express your constraints more generally. Let's look at the first one as an example: mu1 < mu2 < mu3. First, we break it into two constraints: mu1 < mu2 and mu2 < mu3. Then rearrange each:

1. 0 < -mu1 + mu2
2. 0 < -mu2 + mu3

Rewriting as a matrix:

$$\begin{bmatrix} 0 \\ 0 \end{bmatrix} < \begin{bmatrix} -1 & 1 & 0 \\ 0 & -1 & 1 \end{bmatrix} \cdot \begin{bmatrix} mu1 \\ mu2 \\ mu3 \end{bmatrix} < \begin{bmatrix} \infty \\ \infty \end{bmatrix}$$

The above is just for the one constraint - you'll have more rows and variables to express your additional constraints.

You then translate the above into Numpy matrices, and provide it to LinearConstraint.

• First off, thank you so much for this example! It's really helpful! For your example, and looking at LinearConstraint tutorial in Optimization I would write the first linear constraint as : linear_cons=LinearConstraint([[-1,1,0],[0,-1,1]], [0,0], [np.inf,np.inf]. I do the same for the rest. Again, thank you very much. – ivan199415 Feb 22 at 12:35
• I have an additional questions, if you don't mind. 1) I've combined constraints (m1<m2<m3) and m1<=2m3 into one matrix and got [ [-2,1-2],[-1,-1,3] ], is this viable enough? 2) I have no idea how to create a constraint to set w1=w3, specifically using Linear Constraints. I've opted to specify that w1 is within 0.9 to 1.1 of w3. But I still have the problem of setting weights to equal w, w1+w2+w3=1 using linear constraint class. – ivan199415 Feb 22 at 12:57
• I'm not sure if  [ [-2,1-2],[-1,-1,3] ] is right, I haven't checked mathematically, but it doesn't look right. Can you explain how you got that? (maybe edit your original question). w1=w3 => 0 <= w1-w3 <=0 so you can set the lb and ub to be 0. – Cam.Davidson.Pilon Feb 22 at 13:18
• Editted, I've had error, so I think it should be fine now. – ivan199415 Feb 22 at 13:51