python中fit()函数,在python中使用scipy的curvefit拟合Boxcar函数的问题

I can't get this boxcar fit working...I get " OptimizeWarning: Covariance of the parameters could not be estimated

category=OptimizeWarning)", and the output coefficients are not improved beyond the starting guess.

import numpy as np

from scipy.optimize import curve_fit

def box(x, *p):

height, center, width = p

return height*(center-width/2 < x)*(x < center+width/2)

x = np.linspace(-5,5)

y = (-2.5

coeff, var_matrix = curve_fit(box, x, y, p0=[1,0,2])

The output coefficients are [ 1.04499699, 0., 2.], not that the third one has not even been changed.

I suspect that this functional form is not amenable to the levenberg-marquardt algorithm used by curve_fit, which is kind of annoying because I like this function. In mathematica it would be trivial to specify a monte carlo optimization, instead.

解决方案

I suspect that this functional form is not amenable to the levenberg-marquardt algorithm used by curve_fit

You are right. Generally, gradient-based optimizations are not well suited for functions with sharp edges. The gradient is estimated by perturbing the function parameters just a little and looking at the change in fitting quality. However, moving an edge just a little results in zero gradient if it does not cross a data point:

A: it is easy to fit the amplitude because a small change in height immediaterly leads to a change in the residuals.

B: it is hard to fit edges because a small change in position does not affect the residuals (unless the change is big enough to make the edge cross a data point).

Using a stochastic method should work better. Scipy has the differential_evolution function, which uses genetic algorithm and is therefore related to monte-carlo methods. However, it is less trivial to use than curve_fit. You need to specify a cost function and ranges for the parameters:

res = differential_evolution(lambda p: np.sum((box(x, *p) - y)**2), # quadratic cost function

[[0, 2], [-5, 5], [0.1, 10]]) # parameter bounds

It's still a one-liner :)

coeff, var_matrix = curve_fit(box, x, y, p0=[1,0,2])

res = differential_evolution(lambda p: np.sum((box(x, *p) - y)**2), [[0, 2], [-5, 5], [0.1, 10]])

plt.step(x, box(x, *coeff), where='mid', label='curve_fit')

plt.step(x, box(x, *res.x), where='mid', label='diff-ev')

plt.plot(x, y, '.')

plt.legend()