3

我正在尝试在 python 中对具有三个变量的已知函数执行最小二乘拟合。对于随机生成的有错误的数据，我能够完成此任务，但我需要拟合的实际数据包括一些数据点，这些数据点是值的上限。该函数将通量描述为波长的函数，但在某些情况下，在给定波长处测量的通量不是有误差的绝对值，而是通量的最大值，实际值低于该值直至零.

有没有办法告诉拟合任务某些数据点是上限？此外，我必须对许多数据集执行此操作，并且每个数据集可能是上限的数据点数量不同，因此能够自动执行此操作将是有益的，但不是必需的。

如果有任何不清楚的地方，我深表歉意，如果需要，我会努力更清楚地解释。

我用来拟合我的数据的代码包含在下面。

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


def f_all(x,p):
    return np.exp(p[0])/((x**(3+p[1]))*((np.exp(14404.5/((x*1000000)*p[2])))-1))

def residual(p,y,x,error):
    err=(y-(f_all(x,p)))/error
    return err


p0=[-30,2.0,35.0]

data=np.genfromtxt("./Data_Files/Object_001")
wavelength=data[:,0]
flux=data[:,1]
errors=data[:,2]

p,cov,infodict,mesg,ier=leastsq(residual, p0, args = (flux, wavelength, errors), full_output=True)

print p
4

1 回答 1

6

Scipy.optimize.leastsq是一种拟合数据的便捷方式，但下面的工作是函数的最小化。Scipy.optimize包含许多最小化函数，其中一些具有处理约束的能力。这里我解释一下fmin_slsqp我知道的，也许其他人也可以；见Scipy.optimize 文档

fmin_slsqp需要一个函数来最小化和参数的初始值。最小化的函数是残差的平方和。对于参数，我首先执行传统的最小平方拟合，并将结果用作约束最小化问题的初始值。然后有几种施加约束的方法（参见文档）；更简单的是f_ieqcons参数：它需要一个函数，该函数返回一个数组，其值必须始终为正（这就是约束）。如果对于所有最大值点，拟合函数低于该点，则该函数返回正值。

import numpy
import scipy.optimize as scimin
import matplotlib.pyplot as mpl

datax=numpy.array([1,2,3,4,5]) # data coordinates
datay=numpy.array([2.95,6.03,11.2,17.7,26.8])
constraintmaxx=numpy.array([0]) # list of maximum constraints
constraintmaxy=numpy.array([1.2])

# least square fit without constraints
def fitfunc(x,p): # model $f(x)=a x^2+c
    a,c=p
    return c+a*x**2
def residuals(p): # array of residuals
    return datay-fitfunc(datax,p)
p0=[1,2] # initial parameters guess
pwithout,cov,infodict,mesg,ier=scimin.leastsq(residuals, p0,full_output=True) #traditionnal least squares fit

# least square fir with constraints
def sum_residuals(p): # the function we want to minimize
    return sum(residuals(p)**2)
def constraints(p): # the constraints: all the values of the returned array will be >=0 at the end
    return constraintmaxy-fitfunc(constraintmaxx,p)
pwith=scimin.fmin_slsqp(sum_residuals,pwithout,f_ieqcons=constraints) # minimization with constraint

# plotting
ax=mpl.figure().add_subplot(1,1,1)
ax.plot(datax,datay,ls="",marker="x",color="blue",mew=2.0,label="Datas")
ax.plot(constraintmaxx,constraintmaxy,ls="",marker="x",color="red",mew=2.0,label="Max points")
morex=numpy.linspace(0,6,100)
ax.plot(morex,fitfunc(morex,pwithout),color="blue",label="Fit without constraints")
ax.plot(morex,fitfunc(morex,pwith),color="red",label="Fit with constraints")
ax.legend(loc=2)
mpl.show()

在这个例子中，我在抛物线上拟合了一个假想的点样本。这是没有和有约束的结果（左侧的红十字）： 拟合结果

我希望这对您的数据样本有用；否则，请发布您的数据文件之一，以便我们可以尝试使用真实数据。我知道我的示例没有处理数据上的误差线，但您可以通过修改残差函数轻松处理它们。

于 2014-01-09T07:57:34.150 回答