1

I'm trying to solve for the parameters in a mixed gaussian through the method of moments in python.

I've got my set up from the wikipage

I've set up the system of equations in a function, but when I use fsolve from scipy, I'm having issues with the output, after each run I get an error and the output is different.



def myFunction(var):
   (mu_1, sigm_1, mu_2, sigm_2, p) = var
   N = 10_000

   data1 = np.random.normal(1,.5,size = int(N/3))
   data2 = np.random.normal(-1,.5,size = int(N*2/3)+1)
   data = np.concatenate((data1,data2))


   F = np.empty((5))
   F[0] = p*mu_1 + (1-p)*mu_2  +  0
   F[1] = p*(mu_1**2 + sigm_1**2) + (1-p) * (mu_2**2 + sigm_2**2)  + np.std(data)**2
   F[2] = p*(mu_1**3 + 3*mu_1*sigm_1**2) + \
          (1-p) * (mu_2**3 + 3*mu_2 * sigm_2**2) + 0
   F[3] = p*(mu_1**4 + 6*mu_1**2 * sigm_1**2 + 3*sigm_1**2) + \
          (1-p) * (mu_2**4 + 6*mu_2**2 * sigm_2**2 + 3*sigm_2**4) + 3*np.std(data)**4
   F[4] = p*(mu_1**5 + 10*mu_1**3 * sigm_1**2 + 15*mu_1 * sigm_1**4) + \
          (1-p) * (mu_2**5 + 10*mu_2**3 * sigm_2**2 + 15*mu_2*sigm_2**4) + 0

   return F

z_guess = [1,1,1,1,1]
z = sp.optimize.fsolve(myFunction,z_guess)
print(z)

and the error I receive is the following

RuntimeWarning: The iteration is not making good progress, as measured by the 
  improvement from the last five Jacobian evaluations.
  warnings.warn(msg, RuntimeWarning)

I'm not quite sure where my issue is, maybe there's something that I'm missing.

Luis
  • 41
  • 1
  • 5
  • Please include imports in your minimal example. Does this answer your question? https://stackoverflow.com/questions/13057022/python-scipy-optimize-using-fsolve-with-multiple-first-guesses – Natan Mar 17 '22 at 18:40
  • 1
    As you can see, this isn't an error, but a warning, so you can probably ignore it and see what happens. The initial guess is also important: "bad" guesses can lead to poor convergence and even completely different results. Also note that gaussian mixtures are [notoriously difficult to estimate](https://arxiv.org/abs/2202.06930) with the method of moments. It's much simpler to use maximum likelihood estimation via the EM algorithm, although it also has a whole lot of shortcomings. – ForceBru Mar 17 '22 at 18:46

0 Answers0