【问题标题】:What's the error of numpy.polyfit?numpy.polyfit 的错误是什么?
【发布时间】:2013-03-21 04:56:18
【问题描述】:

我想使用numpy.polyfit 进行物理计算,因此我需要误差的大小。

【问题讨论】:

    标签: python numpy curve-fitting


    【解决方案1】:

    如果您在调用polyfit 时指定full=True,它将包含额外信息:

    >>> x = np.arange(100)
    >>> y = x**2 + 3*x + 5 + np.random.rand(100)
    >>> np.polyfit(x, y, 2)
    array([ 0.99995888,  3.00221219,  5.56776641])
    >>> np.polyfit(x, y, 2, full=True)
    (array([ 0.99995888,  3.00221219,  5.56776641]), # coefficients
     array([ 7.19260721]), # residuals
     3, # rank
     array([ 11.87708199,   3.5299267 ,   0.52876389]), # singular values
     2.2204460492503131e-14) # conditioning threshold
    

    返回的残差是拟合误差的平方和,不知道是不是你要的:

    >>> np.sum((np.polyval(np.polyfit(x, y, 2), x) - y)**2)
    7.1926072073491056
    

    在 1.7 版中还有一个 cov 关键字,它将返回系数的协方差矩阵,您可以使用它来计算拟合系数本身的不确定性。

    【讨论】:

    • 你知道 np.polyfit 是使用 TLS(总最小二乘,也称为正交最小二乘)还是 OLS(普通最小二乘)来找到最佳拟合?
    • np.polyfit 使用与scipy.linalg.lstsq 类似的np.linalg.lstsq 的普通最小二乘算法。 scipy.odr 实现正交最小二乘算法,或更准确地说是正交距离回归。
    【解决方案2】:

    正如您在documentation 中看到的:

    Returns
    -------
    p : ndarray, shape (M,) or (M, K)
        Polynomial coefficients, highest power first.
        If `y` was 2-D, the coefficients for `k`-th data set are in ``p[:,k]``.
    
    residuals, rank, singular_values, rcond : present only if `full` = True
        Residuals of the least-squares fit, the effective rank of the scaled
        Vandermonde coefficient matrix, its singular values, and the specified
        value of `rcond`. For more details, see `linalg.lstsq`.
    

    这意味着如果您可以进行拟合并获得残差:

     import numpy as np
     x = np.arange(10)
     y = x**2 -3*x + np.random.random(10)
    
     p, res, _, _, _ = numpy.polyfit(x, y, deg, full=True)
    

    然后,p 是您的拟合参数,res 将是残差,如上所述。 _ 是因为您不需要保存最后三个参数,因此您可以将它们保存在您不会使用的变量 _ 中。这是一个约定,不是必需的。


    @Jaime 的回答解释了残差的含义。您可以做的另一件事是将这些平方偏差视为一个函数(其总和为res)。这对于查看不适合的趋势特别有帮助。 res 可能会很大,因为统计噪声或系统拟合不佳,例如:

    x = np.arange(100)
    y = 1000*np.sqrt(x) + x**2 - 10*x + 500*np.random.random(100) - 250
    
    p = np.polyfit(x,y,2) # insufficient degree to include sqrt
    
    yfit = np.polyval(p,x)
    
    figure()
    plot(x,y, label='data')
    plot(x,yfit, label='fit')
    plot(x,yfit-y, label='var')
    

    所以在图中,注意x = 0附近的不合适:

    【讨论】:

      猜你喜欢
      • 2016-11-11
      • 2018-03-05
      • 2015-05-18
      • 1970-01-01
      • 2021-10-26
      • 2012-06-08
      • 2015-03-28
      • 2017-01-26
      相关资源
      最近更新 更多