python - numpy中的多项式拟合问题

Question

我试图实现一个标准的多项式拟合程序，遇到了一个我无法理解的问题。这是我有样本 x 和 y 数据的代码，我使用正规方程拟合这些数据，还使用 ​​numpy 中的 polyfit 函数。

def mytest(N):
    print "\nPolynomial order (N): {}".format(N)
    mVals = [5, 10, 15, 20]
    a_orig = [198.764, 13.5, 0.523]
    for M in mVals:
        x = arange(-M, M+1)
        y = matrix(a_orig[0]+ a_orig[1]*x + a_orig[2]*x**2).T

        # Code implementing the solution from the normal equations
        nArray = arange(N+1)    
        A = matrix([[n**i for i in nArray] for n in x])
        B = (A.T*A).I
        a_myfit = B*(A.T*y)

        # numpy's polyfit    
        a_polyfit = polyfit(x, y, N)

        print "M: {}".format(M)
        print ["{0:0.3f}".format(i) for i in a_orig] 
        print ["{0:0.3f}".format(i) for i in array(a_myfit)[:,0]]
        print ["{0:0.3f}".format(i) for i in list(array(a_polyfit)[:,0])[::-1]]

mytest(N=5)
mytest(N=6)

这是输出

Polynomial order (N): 5
M: 5
['198.764', '13.500', '0.523']
['198.764', '13.500', '0.523', '0.000', '0.000', '-0.000']
['198.764', '13.500', '0.523', '-0.000', '0.000', '0.000']
M: 10
['198.764', '13.500', '0.523']
['198.764', '13.500', '0.523', '-0.000', '0.000', '-0.000']
['198.764', '13.500', '0.523', '0.000', '0.000', '0.000']
M: 15
['198.764', '13.500', '0.523']
['198.764', '13.500', '0.523', '0.000', '0.000', '-0.000']
['198.764', '13.500', '0.523', '-0.000', '-0.000', '0.000']
M: 20
['198.764', '13.500', '0.523']
['198.764', '13.500', '0.523', '0.000', '-0.000', '-0.000']
['198.764', '13.500', '0.523', '0.000', '0.000', '0.000']

Polynomial order (N): 6
M: 5
['198.764', '13.500', '0.523']
['198.764', '13.500', '0.523', '0.000', '0.000', '-0.000', '-0.000']
['198.764', '13.500', '0.523', '-0.000', '0.000', '0.000', '-0.000']
M: 10
['198.764', '13.500', '0.523']
['198.764', '13.500', '0.523', '-0.000', '-0.000', '-0.000', '-0.000']
['198.764', '13.500', '0.523', '0.000', '0.000', '-0.000', '-0.000']
M: 15
['198.764', '13.500', '0.523']
['294.451', '13.500', '-0.061', '0.000', '-0.001', '-0.000', '-0.000']
['198.764', '13.500', '0.523', '0.000', '0.000', '0.000', '-0.000']
M: 20
['198.764', '13.500', '0.523']
['369.135', '13.500', '-0.046', '0.000', '-0.000', '-0.000', '-0.000']
['198.764', '13.500', '0.523', '0.000', '0.000', '-0.000', '-0.000']

多项式拟合的值给出了 N > 5 和 M > 13 的错误值。

我哪里错了？ polyfit 实现有何不同？

Answer 1

A 是 dtype int32。

可表示的最大值为:

In [10]: np.iinfo(np.dtype('int32')).max
Out[10]: 2147483647

当指数非常大时，n**i 可能会大于 2147483647。此时您会得到不正确的结果。

NumPy does not check for arithmetic overflows (在对数组执行操作时)因为这会影响性能。您可以选择避免算术溢出的数据类型。

要解决这个问题，请声明一个可以表示更大数字的不同数据类型:

A = np.matrix([[n**i for i in nArray] for n in x], dtype='float64')

请注意，问题仍然存在，它只会在 N 的值大得多时发生。