algorithm - 如何计算最适合一组给定数据的二次方程

Question

我有一个包含 20 个实数的向量 X 和一个包含 20 个实数的向量 Y。

我想将它们建模为

y = ax^2+bx + c

如何找到 'a' 、 'b' 和 'c' 的值 和最佳拟合二次方程。

给定值

X = (x1,x2,...,x20)
Y = (y1,y2,...,y20)

我需要一个公式或过程来找到以下值

a = ???
b = ???
c = ???

提前致谢。

Answer 1

@Bartoss 说的都对，+1。我想我只是在这里添加一个实际的实现，没有 QR 分解。您想要评估 a、b、c 的值，以使测量数据和拟合数据之间的距离最小。您可以选择作为衡量标准

sum(ax^2+bx + c -y)^2)

其中对向量 x,y 的元素求和。

那么，最小值意味着数量对 a、b、c 中每一个的导数都为零:

d (sum(ax^2+bx + c -y)^2) /da =0
d (sum(ax^2+bx + c -y)^2) /db =0
d (sum(ax^2+bx + c -y)^2) /dc =0

这些方程是

2(sum(ax^2+bx + c -y)*x^2)=0
2(sum(ax^2+bx + c -y)*x)  =0
2(sum(ax^2+bx + c -y))    =0

除以2，上面可以改写为

a*sum(x^4) +b*sum(x^3) + c*sum(x^2) =sum(y*x^2)
a*sum(x^3) +b*sum(x^2) + c*sum(x)   =sum(y*x)
a*sum(x^2) +b*sum(x)   + c*N        =sum(y)

在你的情况下，N=20。下面是一个简单的 Python 代码，展示了如何执行此操作。

from numpy import random, array
from scipy.linalg import solve
import matplotlib.pylab as plt
a, b, c = 6., 3., 4.
N = 20
x = random.rand((N))
y = a * x ** 2 + b * x + c
y += random.rand((20)) #add a bit of noise to make things more realistic

x4 = (x ** 4).sum()
x3 = (x ** 3).sum()
x2 = (x ** 2).sum()
M = array([[x4, x3, x2], [x3, x2, x.sum()], [x2, x.sum(), N]])
K = array([(y * x ** 2).sum(), (y * x).sum(), y.sum()])
A, B, C = solve(M, K)

print 'exact values     ', a, b, c
print 'calculated values', A, B, C

fig, ax = plt.subplots()
ax.plot(x, y, 'b.', label='data')
ax.plot(x, A * x ** 2 + B * x + C, 'r.', label='estimate')
ax.legend()
plt.show()

实现解决方案的一种更快的方法是使用非线性最小二乘算法。这样写起来会更快，但运行起来不会更快。使用 scipy 提供的，

from scipy.optimize import leastsq
def f(arg):
    a,b,c=arg
    return a*x**2+b*x+c-y

(A,B,C),_=leastsq(f,[1,1,1])#you must provide a first guess to start with in this case.