Lorentzian scipy.optimize.leastsq 拟合数据失败

Question

自从我上了Python的课程后，我就想用它来处理我的数据。虽然我尝试了很久，但还是不知道为什么这段代码不管用。

我想做的事情

我想从一个名为“Test”的子文件夹中逐个读取数据文件，对数据进行一些处理，然后用洛伦兹函数来拟合这些数据。

问题描述

当我运行下面的代码时，它什么都没有拟合，只是在调用了4次函数后返回了我最初的参数。我尝试过缩放数据，调整过ftol和maxfev这两个参数，反复查看Python文档，但没有任何改善。我还尝试将列表明确地转换为numpy.arrays，以及参考了一个问题的解决方案scipy.optimize.leastsq返回最佳猜测参数而不是新的最佳拟合，使用了x = x.astype(np.float64)，但也没有效果。奇怪的是，对于少数几个特定的数据文件，这段代码在某个时候是有效的，但对于大多数文件来说，它从来没有成功过。可以肯定的是，这些数据是可以拟合的，因为在Origin中使用Levenberg-Marquard拟合程序可以得到相当不错的结果。

有人能告诉我哪里出错了吗，或者给我一些其他的建议吗...?

import numpy,math,scipy,pylab
from scipy.optimize import leastsq
import glob,os
for files in glob.glob("*.txt"):
    x=[]
    y=[]
    z=[]
    f = open(files, 'r')
    raw=f.readlines()
    f.close()
    del raw[0:8]       #delete Header
    for columns in ( raw2.strip().split() for raw2 in raw ):  #data columns
        x.append(float(columns[0]))
        y.append(float(columns[1]))
        z.append(10**(float(columns[1])*0.1)) #transform data for the fit
    def lorentz(p,x):
        return (1/(1+(x/p[0] - 1)**4*p[1]**2))*p[2]
    def errorfunc(p,x,z):
        return lorentz(p,x)-z

    p0=[3.,10000.,0.001]

    Params,cov_x,infodict,mesg,ier = leastsq(errorfunc,p0,args=(x,z),full_output=True)
    print Params
    print ier

Answer 1

在没有看到你的数据之前，很难判断问题出在哪里。我生成了一些随机噪声，并用你的代码对它进行了拟合。结果一切正常。这个算法不允许设置参数的边界，所以如果你的 p0 值接近零，可能会遇到问题。我做了以下操作：

import numpy as np
from scipy.optimize import leastsq
import matplotlib.pyplot as plt

def lorentz(p,x):
    return p[2] / (1.0 + (x / p[0] - 1.0)**4 * p[1]**2)

def errorfunc(p,x,z):
        return lorentz(p,x)-z

p = np.array([0.5, 0.25, 1.0], dtype=np.double)
x = np.linspace(-1.5, 2.5, num=30, endpoint=True)
noise = np.random.randn(30) * 0.05
z = lorentz(p,x)
noisyz = z + noise

p0 = np.array([-2.0, -4.0, 6.8], dtype=np.double) #Initial guess
solp, ier = leastsq(errorfunc, 
                    p0, 
                    args=(x,noisyz),
                    Dfun=None,
                    full_output=False,
                    ftol=1e-9,
                    xtol=1e-9,
                    maxfev=100000,
                    epsfcn=1e-10,
                    factor=0.1)

plt.plot(x, z, 'k-', linewidth=1.5, alpha=0.6, label='Theoretical')
plt.scatter(x, noisyz, c='r', marker='+', color='r', label='Measured Data')
plt.plot(x, lorentz(solp,x), 'g--', linewidth=2, label='leastsq fit')
plt.xlim((-1.5, 2.5))
plt.ylim((0.0, 1.2))
plt.grid(which='major')
plt.legend(loc=8)
plt.show()

这得到了一个解：
solp = array([ 0.51779002, 0.26727697, 1.02946179])
这个解接近理论值：
np.array([0.5, 0.25, 1.0])