如何在python中使用仅为上限的数据进行最小二乘拟合？

Question

我正在尝试在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

Answer 1

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

fmin_slsqp需要最小化函数和参数的初始值。最小化函数是残差平方和。对于参数，我首先执行传统的leatsq拟合，并将结果用作约束最小化问题的初始值。然后有几种方法来施加约束（参见doc）；更简单的是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()

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

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