python - Lorentzian scipy.optimize.leastsq 适合数据失败

Question

自从我参加了 Python 讲座后，我想用它来拟合我的数据。虽然我现在已经尝试了一段时间，但我仍然不知道为什么这不起作用。

我想做什么

从子文件夹(此处称为“测试”)中取出一个又一个数据文件，稍微转换数据并用洛伦兹函数对其进行拟合。

问题描述

当我运行下面发布的代码时，它不适合任何东西，在 4 次函数调用后只返回我的初始参数。在一遍又一遍地检查 python 文档后，我尝试缩放数据，尝试使用 ftol 和 maxfev，但没有任何改进。我还尝试将列表显式更改为 numpy.arrays，以及对问题 scipy.optimize.leastsq returns best guess parameters not new best fit 的解决方案, x = x.astype(np.float64)。没提升。奇怪的是，对于少数选定的数据文件，同样的代码在某些时候起作用，但对于大多数人来说它从来没有起作用。它绝对可以拟合，因为 Levenberg-Marquard 拟合例程在 Origin 中给出了相当好的结果。

有人能告诉我哪里出了问题或指出替代方案...？

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])