8

自从我上过关于 Python 的讲座以来，我想用它来拟合我的数据。尽管我已经尝试了一段时间，但我仍然不知道为什么这不起作用。

我想做什么

从子文件夹（此处称为“测试”）中获取一个接一个的数据文件，稍微转换数据并使其与洛伦兹函数相匹配。

问题描述

当我运行下面发布的代码时，它不适合任何东西，只是在 4 次函数调用后返回我的初始参数。我尝试缩放数据，ftol反复maxfev检查 python 文档，但没有任何改进。我还尝试将列表更改为numpy.arrays明确，以及为问题scipy.optimize.leastsq 返回最佳猜测参数而不是新的最佳拟合给出的解决方案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
4

1 回答 1

7

如果不查看您的数据，就很难判断出了什么问题。我产生了一些随机噪音，并使用您的代码对其进行了拟合。一切正常。该算法不允许参数边界，因此如果您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]) 在此处输入图像描述

于 2013-08-28T18:25:58.200 回答