Why does scipy.optimize.curve_fit produce parameters which are barely different from the guess?

Question

I've been trying to fit some histogram data with scipy.optimize.curve_fit, but so far I haven't once been able to produce fit parameters that differ significantly from my guess parameters.

I wouldn't be terribly surprised to find that the more arcane parameters in my fit get stuck in local minima, but even linear coefficients won't move from my initial guesses!

If you've seen anything like this before, I'd love some advice. Do least-squared minimization routines just not work for certain classes of functions?

I try this,

import numpy as np
from matplotlib.pyplot import *
from scipy.optimize import curve_fit

def grating_hist(x,frac,xmax,x0):
  #  model data to be turned into a histogram
  dx = x[1]-x[0]
  z = np.linspace(0,1,20000,endpoint=True)
  grating = np.cos(frac*np.pi*z)
  norm_grating = xmax*(grating-grating[-1])/(1-grating[-1])+x0
  # produce the histogram
  bin_edges = np.append(x,x[-1]+x[1]-x[0])
  hist,bin_edges = np.histogram(norm_grating,bins=bin_edges)
  return hist

x = np.linspace(0,5,512)
p_data = [0.7,1.1,0.8]
pct = grating_hist(x,*p_data)
p_guess = [1,1,1]
p_fit,pcov = curve_fit(grating_hist,x,pct,p0=p_guess)

plot(x,pct,label='Data')
plot(x,grating_hist(x,*p_fit),label='Fit')
legend()

show()

print 'Data Parameters:', p_data
print 'Guess Parameters:', p_guess
print 'Fit Parameters:', p_fit
print 'Covariance:',pcov

and I see this: http://i.stack.imgur.com/GwXzJ.png (I'm new here, so I can't post images)

Data Parameters: [0.7, 1.1, 0.8]
Guess Parameters: [1, 1, 1]
Fit Parameters: [ 0.97600854  0.99458336  1.00366634]
Covariance: [[  3.50047574e-06  -5.34574971e-07   2.99306123e-07]
 [ -5.34574971e-07   9.78688795e-07  -6.94780671e-07]
 [  2.99306123e-07  -6.94780671e-07   7.17068753e-07]]

Whaaa? I'm pretty sure this isn't a local minimum for variations in xmax and x0, and it's a long way from the global minimum best fit. The fit parameters still don't change, even with better guesses. Different choices for curve functions (e.g. the sum of two normal distributions) do produce new parameters for the same data, so I know it's not the data itself. I also tried the same thing with scipy.optimize.leastsq itself just in case, but no dice; the parameters still don't move. If you have any thoughts on this, I'd love to hear them!

Answer 1

The problem you're facing is actually not due to curve_fit (or leastsq). It is due to the landscape of the objective of your optimisation problem. In your case the objective is the sum of residuals' squares, which you are trying to minimise. Now, if you look closely at your objective in a close surrounding of your initial conditions, for example using the code below, which only focuses on the first parameter:

p_ind = 0
eps = 1e-6
n_points = 100

frac_surroundings = np.linspace(p_guess[p_ind] - eps, p_guess[p_ind] + eps, n_points)

obj = []
temp_guess = p_guess.copy()
for p in frac_surroundings:
    temp_guess[0] = p
    obj.append(((grating_hist(x, *p_data) - grating_hist(x, *temp_guess))**2.0).sum())

py.plot(frac_surroundings, obj)
py.show()

you will notice that the landscape is piecewise constant (you can easily check that the situation is the same for other parameters. The problem with that is that these pieces are of the order of 10^-6, whereas the initial step of the fitting procedure is somewhere around 10^-8, hence the procedure ends quickly concluding that you cannot improve from the given initial condition. You could try to fix it by changing epsfcn parameter in curve_fit, but you would quickly notice that the landscape, on top of being piecewise constant, is also very "rugged". In other words, curve_fit is simply not well suited for such a problem, which is simply difficult for gradient based methods, as it is highly non-convex. Probably, some stochastic optimisation methods could do a better job. That is, however, a different question/problem.

Answer 2

I think it is a local minimum, or the algorith fails for a non trivial reason. It is far easier to fit the data to the input, instead of fitting the statistical description of the data to the statistical description of the input.

Here's a modified version of the code doing so:

z = np.linspace(0,1,20000,endpoint=True)

def grating_hist_indicator(x,frac,xmax,x0):
  #  model data to be turned into a histogram
  dx = x[1]-x[0]
  grating = np.cos(frac*np.pi*z)
  norm_grating = xmax*(grating-grating[-1])/(1-grating[-1])+x0
  return norm_grating

x = np.linspace(0,5,512)
p_data = [0.7,1.1,0.8]
pct = grating_hist(x,*p_data)

pct_indicator = grating_hist_indicator(x,*p_data)
p_guess = [1,1,1]
p_fit,pcov = curve_fit(grating_hist_indicator,x,pct_indicator,p0=p_guess)

plot(x,pct,label='Data')
plot(x,grating_hist(x,*p_fit),label='Fit')
legend()
show()