python - 如何使用仅为上限的数据在 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 需要一个最小化函数和参数的初始值。要最小化的函数是残差平方和。对于参数，我首先执行传统的 leastsq 拟合，并将结果用作约束最小化问题的初始值。然后有几种方法可以施加约束(参见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()

在这个例子中，我在抛物线上拟合了一个假想的点样本。这是没有和有约束的结果(左边的红叉):

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