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.
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]
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:
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.