python - 将结果与二次回归混淆

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

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

NumPy 之所以没有这个问题，是因为他们用最小二乘法求解系统，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%

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