将结果与二次回归混淆

Question

所以，我正在尝试用二次回归拟合一些 x,y 数据对，可以在 http://polynomialregression.drque.net/math.html 找到示例公式。以下是我的代码，它使用该显式公式并使用 numpy 内置函数进行回归，

import numpy as np 
x = [6.230825,6.248279,6.265732]
y = [0.312949,0.309886,0.306639472]
toCheck = x[2]


def evaluateValue(coeff,x):
    c,b,a = coeff
    val = np.around( a+b*x+c*x**2,9)
    act = 0.306639472
    error=  np.abs(act-val)*100/act
    print "Value = {:.9f} Error = {:.2f}%".format(val,error)



###### USing numpy######################
coeff = np.polyfit(x,y,2)
evaluateValue(coeff, toCheck)



################# Using explicit formula
def determinant(a,b,c,d,e,f,g,h,i):
    # the matrix is [[a,b,c],[d,e,f],[g,h,i]]
    return a*(e*i - f*h) - b*(d*i - g*f) + c*(d*h - e*g)

a = b = c = d = e = m = n = p = 0
a = len(x)
for i,j in zip(x,y):
        b += i
        c += i**2
        d += i**3
        e += i**4
        m += j
        n += j*i
        p += j*i**2
det = determinant(a,b,c,b,c,d,c,d,e)
c0 = determinant(m,b,c,n,c,d,p,d,e)/det
c1 = determinant(a,m,c,b,n,d,c,p,e)/det
c2 = determinant(a,b,m,b,c,n,c,d,p)/det

evaluateValue([c2,c1,c0], toCheck)




######Using another explicit alternative
def determinantAlt(a,b,c,d,e,f,g,h,i):
    return a*e*i - a*f*h - b*d*i +b*g*f + c*d*h - c*e*g # <- barckets removed

a = b = c = d = e = m = n = p = 0
a = len(x)
for i,j in zip(x,y):
        b += i
        c += i**2
        d += i**3
        e += i**4
        m += j
        n += j*i
        p += j*i**2
det = determinantAlt(a,b,c,b,c,d,c,d,e)
c0 = determinantAlt(m,b,c,n,c,d,p,d,e)/det
c1 = determinantAlt(a,m,c,b,n,d,c,p,e)/det
c2 = determinantAlt(a,b,m,b,c,n,c,d,p)/det

evaluateValue([c2,c1,c0], toCheck)

此代码给出此输出

Value = 0.306639472 Error = 0.00%
Value = 0.308333580 Error = 0.55%
Value = 0.585786477 Error = 91.03%

作为，你可以看到这些是彼此不同的，第三个是完全错误的。现在我的问题是：
1. 为什么显式公式给出的结果有点错误，如何改进？
2. numpy 如何给出如此准确的结果？
3.第三种情况只是打开括号，结果怎么变化这么大？

Answer 1

所以这里发生的一些事情不幸地困扰着你做事的方式。看看这段代码：

for i,j in zip(x,y):
        b += i
        c += i**2
        d += i**3
        e += i**4
        m += j
        n += j*i
        p += j*i**2

您构建的特征使得 x 值不仅是平方的，而且是立方和四次方的。

如果在将这些值放入 3 x 3 矩阵以求解之前打印出每个值：

In [35]: a = b = c = d = e = m = n = p = 0
    ...: a = len(x)
    ...: for i,j in zip(xx,y):
    ...:         b += i
    ...:         c += i**2
    ...:         d += i**3
    ...:         e += i**4
    ...:         m += j
    ...:         n += j*i
    ...:         p += j*i**2
    ...: print(a, b, c, d, e, m, n, p)
    ...:
    ...:
3 18.744836 117.12356813829001 731.8283056811686 4572.738547313946 0.9294744720000001 5.807505391292503 36.28641270376207

在处理 floating-point 算术时，尤其是对于小值，运算顺序确实很重要。这里发生的事情是侥幸，计算出的小值和大值的混合导致了一个非常小的值。因此，当您使用分解形式和扩展形式计算行列式时，请注意您如何获得略有不同的结果，但也要查看值的精度：

In [36]: det = determinant(a,b,c,b,c,d,c,d,e)

In [37]: det
Out[37]: 1.0913403514223319e-10

In [38]: det = determinantAlt(a,b,c,b,c,d,c,d,e)

In [39]: det
Out[39]: 2.3283064365386963e-10

行列式在10^-10量级！存在差异的原因是因为使用 floating-point 算术，理论上两种行列式方法应该产生相同的结果，但不幸的是，实际上它们给出的结果略有不同，这是由于称为 错误传播 。因为可以表示 floating-point 数字的位数是有限的，所以运算顺序会改变错误传播的方式，因此即使您删除了括号并且公式基本匹配，运算顺序也会得到结果现在不同了。对于经常处理 floating-point 算术的任何软件开发人员来说，这篇文章都是必读的：What Every Computer Scientist Should Know About Floating-Point Arithmetic.

因此，当您尝试使用 Cramer 规则求解系统时，不可避免地要除以代码中的主要行列式，即使变化大约为 10^-10，这两种方法之间的变化可以忽略不计，但是您会得到非常不同的结果，因为您在求解系数时除以这个数字。

NumPy 没有这个问题的原因是他们通过 least-squares 和 pseudo-inverse 解决了系统，而不是使用 Cramer 规则。我不建议使用 Cramer 规则来查找回归系数，这主要是由于经验，而且还有更可靠的方法。

然而，要解决您的特定问题，最好对数据进行 标准化 以便动态范围现在以 0 为中心。因此，您用来构建系数矩阵的特征更明智，因此计算过程更容易处理数据。在您的情况下，像用 x 值的平均值减去数据这样简单的事情应该可行。因此，如果您有要预测的新数据点，您必须在进行预测之前首先减去 x 数据的平均值。

因此，在您的代码开头，对该数据执行均值减法和回归。我已经向您展示了我在哪里修改了上面给出的源代码：

import numpy as np 
x = [6.230825,6.248279,6.265732]
y = [0.312949,0.309886,0.306639472]

# Calculate mean
me = sum(x) / len(x)
# Make new dataset that is mean subtracted
xx = [pt - me for pt in x]

#toCheck = x[2]

# Data point to check is now mean subtracted
toCheck = x[2] - me



def evaluateValue(coeff,x):
    c,b,a = coeff
    val = np.around( a+b*x+c*x**2,9)
    act = 0.306639472
    error=  np.abs(act-val)*100/act
    print("Value = {:.9f} Error = {:.2f}%".format(val,error))



###### USing numpy######################
coeff = np.polyfit(xx,y,2) # Change
evaluateValue(coeff, toCheck)



################# Using explicit formula
def determinant(a,b,c,d,e,f,g,h,i):
    # the matrix is [[a,b,c],[d,e,f],[g,h,i]]
    return a*(e*i - f*h) - b*(d*i - g*f) + c*(d*h - e*g)

a = b = c = d = e = m = n = p = 0
a = len(x)
for i,j in zip(xx,y): # Change
        b += i
        c += i**2
        d += i**3
        e += i**4
        m += j
        n += j*i
        p += j*i**2
det = determinant(a,b,c,b,c,d,c,d,e)
c0 = determinant(m,b,c,n,c,d,p,d,e)/det
c1 = determinant(a,m,c,b,n,d,c,p,e)/det
c2 = determinant(a,b,m,b,c,n,c,d,p)/det

evaluateValue([c2,c1,c0], toCheck)




######Using another explicit alternative
def determinantAlt(a,b,c,d,e,f,g,h,i):
    return a*e*i - a*f*h - b*d*i +b*g*f + c*d*h - c*e*g # <- barckets removed

a = b = c = d = e = m = n = p = 0
a = len(x)
for i,j in zip(xx,y): # Change
        b += i
        c += i**2
        d += i**3
        e += i**4
        m += j
        n += j*i
        p += j*i**2
det = determinantAlt(a,b,c,b,c,d,c,d,e)
c0 = determinantAlt(m,b,c,n,c,d,p,d,e)/det
c1 = determinantAlt(a,m,c,b,n,d,c,p,e)/det
c2 = determinantAlt(a,b,m,b,c,n,c,d,p)/det
evaluateValue([c2,c1,c0], toCheck)

当我运行这个时，我们现在得到：

In [41]: run interp_test
Value = 0.306639472 Error = 0.00%
Value = 0.306639472 Error = 0.00%
Value = 0.306639472 Error = 0.00%

作为您的最终读物，这是我在他们的问题中解决的其他人遇到的类似问题：。总结就是我建议他们不要使用 Cramer 规则，而是使用 least-squares 到 pseudo-inverse。我向他们展示了如何在不使用 numpy.polyfit 的情况下获得完全相同的结果。此外，使用 least-squares 概括了如果你有超过 3 个点，你仍然可以通过你的点拟合二次方，以便模型具有尽可能小的误差。

将结果与二次回归混淆

Confusing result with quadratic regression

python

regression

numpy

curve-fitting

smoothing