与二次回归混淆的结果

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

在处理浮点运算时，尤其是对于小值，运算顺序很重要。这里发生的情况是，由于侥幸，已计算的小值和大值的混合导致一个非常小的值。因此，当您使用因式分解和展开式计算行列式时，请注意您得到的结果略有不同，但也要注意值的精度：

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

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

NumPy 没有这个问题的原因是他们通过最小二乘法和pseudo-inverse 解决系统问题，而不是使用克莱默规则。我不建议使用 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%

作为您的最后阅读，这是我在他们的问题中解决的其他人遇到的类似问题：Fitting a quadratic function in python without numpy polyfit。总结是我建议他们不要使用克莱默规则并通过伪逆使用最小二乘法。我向他们展示了如何在不使用numpy.polyfit 的情况下获得完全相同的结果。此外，使用最小二乘法可以概括如果您有 3 个以上的点，您仍然可以通过您的点拟合二次方，从而使模型具有尽可能小的误差。