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

Question

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

我想做什么

从子文件夹中获取一个接一个的数据文件(此处称为'Test'),稍微转换数据并使用Lorentzian函数拟合.

问题描述

当我运行下面发布的代码时,它不适合任何东西,只是在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

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