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

Question

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):
 from autograd.scipy.stats import norm
 from autograd import numpy
 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(
value_and_grad(negative_log_likelihood), # see autograd docs.
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.

Answer 1

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.