利用pyhton拟和区域边界

本文翻译自python - 将曲线拟合到散点图的边界 - IT工具网 (coder.work),如有侵权,请联系删除.

一.直接调用

只需修改参数表达式model,参数初始预测值guesses,再根据需要选择拟和上边缘或者下边缘

# 选取需要拟和的列
X = df[['横坐标列']].values
y = df['纵坐标列'].values

#绘制原始数据
plt.figure(dpi=100,figsize=(8,4))
plt.plot(X.T[0] , y, '.', alpha=0.2, label = '原始数据') 

#定义拟和函数
model = lambda x, a, b: a*x+b  #参数表达式
guesses =[200, -150] #参数初始预测值

#定义残差调整函数
def get_flipped(y_data, y_model):
    flipped = y_model - y_data
    flipped[flipped < 0] = 0#上边缘拟和
#     flipped[flipped > 0] = 0#下边缘拟和
    return flipped

#定义残差函数
def flipped_resid(pars, x, y):
    y_model = model(x, *pars)
    flipped = get_flipped(y, y_model)
    resid = np.square(y + flipped - y_model)
    return np.nan_to_num(resid)

#调用leastsq拟和
fit_pars, flag = leastsq(func = flipped_resid, x0 = guesses,
                         args = (X.T[0], y))
y_fit = model(X.T[0], *fit_pars)
plt.plot(X.T[0], y_fit, 'r--', zorder = 0.9)

二.详细探讨

1.正常最小二乘残差函数

def error(pars, x, y):
    y_model = model(x, *pars)#参数代入计算拟和纵坐标
    resid = y - y_model#计算残差
    return np.nan_to_num(resid)#替换空值和无穷值

2.残差调整函数

最小二乘的目的是求残差函数y + flipped - y_model等于0时,对应的参数值.

如果该表达式不为0就继续迭代.

def get_flipped(y_data, y_model):#传入原始纵坐标和拟和纵坐标
    flipped = y_model - y_data#计算假残差
    flipped[flipped < 0] = 0
#如果flipped>=0,即拟和纵坐标大于等于原始纵坐标,也即拟和到达上边缘,
# 返回flipped = y_model - y_data
# 又定义误差函数y-y_model+flipped=y-y_model+(y_model - y_data)=0
# 即到达上边缘时完成最小化误差函数的任务,最小二乘结束并返回参数

#如果flipped<0,即拟和纵坐标小于原始纵坐标,也即拟和未到达拟和上边缘,
# 返回调整后的flipped 不等于 y_model - y_data
# 又定义误差函数y-y_model+flipped不等于0,继续迭代

#综上,最小二乘逐次增大拟和纵坐标,直到到达上边缘
    return flipped

3.残差函数

def flipped_resid(pars, x, y):#传入原数据,以及每次迭代后的参数
    y_model = model(x, *pars)#参数代入计算拟和纵坐标
    flipped = get_flipped(y, y_model)#将原始纵坐标和拟和纵坐标传入,返回调整后flipped
    resid = np.square(y-y_model+flipped)#计算残差,平方可以理解为放大残差进行精确的拟和
    #print pars, resid.sum()# 取消注释可以打印每次迭代的参数和残差平方和
    return np.nan_to_num(resid)#替换空值和无穷值

4.调用leastsq拟和

model = lambda x, a, b: a*x+b  #这里拟和一条直线,根据实际需要更改表达式
guesses =[200, -150]#给出参数a,b初始预测值
#调用leastsq
# func:误差函数
# x0:参数a,b的初始预测值 
# args:拟和点的横纵坐标,np数组形式 
# fit_pars:保存参数a,b最终拟和结果的列表
# flag:可忽略
fit_pars, flag = leastsq(func = flipped_resid, x0 = guesses,
                         args = (X.T[0], y))
y_fit = model(X.T[0], *fit_pars)#最小二乘拟和参数代入,计算拟和纵坐标
fit_pars=np.around(fit_pars,2)#保留2位小数
text_above=('y=%.3f*x+(%.3f)'%(fit_pars[0],fit_pars[1]))#标签
plt.plot(X.T[0], y_fit, 'r--', zorder = 0.9, label = text_above)