编程语言
首页 > 编程语言> > 如何使用仅仅是上限的数据在python中进行最小二乘拟合?

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

作者:互联网

我试图在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

解决方法:

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

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

在这个例子中,我在抛物线上拟合了一个虚构的点样本.这是结果,没有约束(左边是红叉):

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

标签:least-squares,model-fitting,python
来源: https://codeday.me/bug/20190831/1773681.html