如何在 Python 中实现约束线性拟合？

Question

我正在尝试将线性模型拟合到一组数据，约束是所有残差（模型 - 数据）都是正的 - 换句话说，模型应该是“最佳高估”。没有这个约束，线性模型可以很容易地用 numpy 的 polyfit 找到，如下所示。

import numpy as np
import matplotlib.pyplot as plt

x = [-4.12179107e-01, -1.40664082e-01, -5.52301563e-06,  1.82898473e-01]
y = [-4.14846251, -3.31607886, -3.57827245, -5.09914559]

plt.scatter(x,y)
coeff = np.polyfit(x,y,1)
plt.plot(x,np.polyval(coeff,x),c='r',label='numpy-polyval')
plt.plot(x,np.polyval([-2,-3.6],x),c='g',label='desired-fit') #a rough guess of the desired result
plt.legend()

example1

有没有一种有效的方法来实现这种约束的线性拟合？

Answer 1

是的，最合适的是通过顶部两点的线。我做了一个 argsort 来找到顶部的 Y，计算斜率和 y-intercept，然后我们开始：

import numpy as np
import matplotlib.pyplot as plt

x = [-4.12179107e-01, -1.40664082e-01, -5.52301563e-06,  1.82898473e-01]
y = [-4.14846251, -3.31607886, -3.57827245, -5.09914559]

plt.scatter(x,y)
coeff = np.polyfit(x,y,1)
model1 = np.polyval(coeff,x)
model1 += (y-model1).max()
print(model1)
print(sum((y-model1)**2))

z = np.argsort(y)
pt0 = (x[z[-1]],y[z[-1]])
pt1 = (x[z[-2]],y[z[-2]])
m = (pt1[1]-pt0[1])/(pt1[0]-pt0[0])
b = pt0[1]-m*pt0[0]

model2 = np.polyval([m,b],x)
print(model2)
print(sum((y-model2)**2))

plt.plot(x,model1,c='r',label='numpy-polyval')
plt.plot(x,model2,c='g',label='generated')
plt.legend()
plt.show()

输出：

Answer 2

这是一个二次规划问题。有几个库（CVXOPT、quadprog 等）可用于解决它。下面是一个使用 quadprog 的例子：

import numpy as np
import matplotlib.pyplot as plt
import quadprog

x = [-4.12179107e-01, -1.40664082e-01, -5.52301563e-06, 1.82898473e-01]
y = [-4.14846251, -3.31607886, -3.57827245, -5.09914559]

A = np.c_[x, np.ones(len(x))]
y = np.array(y)
G = A.T @ A
a = A.T @ y
C = A.T
b = y
coeffs = quadprog.solve_qp(G, a, C, b)[0]

plt.scatter(x, y)
plt.plot(x, np.polyval(coeffs, x), c='r')
plt.show()

这给出：

参见例如this post 获取更多信息。特别是描述了如何将线性回归问题设置为二次规划问题。

附带说明一下，最佳线总是会通过一个数据点，但不需要通过两个这样的点。例如，取 x = [-1., 0., 1.] 和 y = [1., 2., 1.].

如何在 Python 中实现约束线性拟合？

How to implement a constrained linear fit in Python?

python

numpy

constraints

curve-fitting