Python科学计算——任意波形拟合

Posted 2023-04-24

参考技术A

任意波形的生成 (geneartion of arbitrary waveform) 在商业，军事等领域都有着重要的应用，诸如空间光通信 (free-space optics communication)， 高速信号处理 (high-speed signal processing)，雷达 (radar) 等。在任意波形生成后， 如何评估生成的任意波形 成为另外一个重要的话题。

假设有一组实验数据，已知他们之间的函数关系：y=f(x)，通过这些信息，需要确定函数中的一些参数项。例如，f 是一个线型函数 f(x)=k*x+b，那么参数 k 和 b 就是需要确定的值。如果这些参数用 p 表示的话，那么就需要找到一组 p 值使得如下公式中的 S 函数最小：

这种算法被称之为 最小二乘拟合 (least-square fitting)。scipy 中的子函数库 optimize 已经提供实现最小二乘拟合算法的函数 leastsq 。下面是 leastsq 函数导入的方式：

scipy.optimize.leastsq 使用方法

在 Python科学计算——Numpy.genfromtxt 一文中，使用 numpy.genfromtxt 对数字示波器采集的三角波数据导入进行了介绍，今天，就以 4GHz三角波 波形的拟合为案例介绍任意波形的拟合方法。

在 Python科学计算——如何构建模型？ 一文中，讨论了如何构建三角波模型。在标准三角波波形的基础上添加了 横向，纵向的平移和伸缩特征参数 ，最后添加了 噪声参数 模拟了三角波幅度参差不齐的随机性特征。但在波形拟合时，并不是所有的特征参数都要纳入考量，例如，噪声参数应是 波形生成系统 的固有特征，正因为它的存在使得产生的波形存在瑕疵，因此，在进行波形拟合并评估时，不应将噪声参数纳入考量，最终模型如下：

在调用 scipy.optimize.leastsq 函数时，需要构建误差函数：

有时候，为了使图片有更好的效果，需要对数据进行一些处理：

leastsq 调用方式如下：

合理的设置 p0 可以减少程序运行时间，因此，可以在运行一次程序后，用拟合后的相应数据对 p0 进行修正。

在对波形进行拟合后，调用 pylab 对拟合前后的数据进行可视化：

均方根误差 (root mean square error) 是一个很好的评判标准，它是观测值与真值偏差的平方和观测次数n比值的平方根，在实际测量中，观测次数n总是有限的，真值只能用最可信赖（最佳）值来代替.方根误差对一组测量中的特大或特小误差反映非常敏感，所以，均方根误差能够很好地反映出测量的精密度。

RMSE 用程序实现如下：

拟合效果，模型参数输出：

leastsq 函数适用于任何波形的拟合，下面就来介绍一些常用的其他波形：

如何使用 Python 和 Numpy 计算 r 平方？

【中文标题】如何使用 Python 和 Numpy 计算 r 平方？

【英文标题】：How do I calculate r-squared using Python and Numpy?

【发布时间】：2010-10-27 23:27:39

【问题描述】：

我正在使用 Python 和 Numpy 计算任意次数的最佳拟合多项式。我传递了一个 x 值、y 值以及我想要拟合的多项式的次数（线性、二次等）的列表。

这很有效，但我还想计算 r（相关系数）和 r-squared（确定系数）。我将我的结果与 Excel 的最佳拟合趋势线功能以及它计算的 r 平方值进行比较。使用它，我知道我正在为线性最佳拟合（度数等于 1）正确计算 r 平方。但是，我的函数不适用于度数大于 1 的多项式。

Excel 能够做到这一点。如何使用 Numpy 计算高阶多项式的 r 平方？

这是我的功能：

import numpy

# Polynomial Regression
def polyfit(x, y, degree):
    results = 

    coeffs = numpy.polyfit(x, y, degree)
     # Polynomial Coefficients
    results['polynomial'] = coeffs.tolist()

    correlation = numpy.corrcoef(x, y)[0,1]

     # r
    results['correlation'] = correlation
     # r-squared
    results['determination'] = correlation**2

    return results

【问题讨论】：

注意：您只在计算系数时使用度数。

tydok 是正确的。您正在计算 x 和 y 的相关性以及 y=p_0 + p_1 * x 的 r 平方。请参阅下面的我的答案以获取一些应该工作的代码。如果你不介意我问，你的最终目标是什么？你在做模型选择（选择使用什么程度）？还是别的什么？

@leif -- 请求归结为“像 Excel 那样做”。我从这些答案中感觉到，在使用非线性最佳拟合曲线时，用户可能对 r 平方值的解读过多。尽管如此，我不是数学向导，这是要求的功能。

附带问题：pandas corr() 函数不返回 r^"2 pearson 系数吗？

【参考方案1】：

R-squared 是一种仅适用于线性回归的统计量。

本质上，它衡量的是数据中有多少变化可以用线性回归来解释。

因此，您计算“总平方和”，即每个结果变量与其均值的总平方偏差。 . .

其中 y_bar 是 y 的平均值。

然后，您计算“回归平方和”，即您的 FITTED 值与平均值相差多少

并找出这两者的比率。

现在，多项式拟合所需要做的就是插入该模型中的 y_hat，但称其为 r-squared 并不准确。

Here 是我发现的一个链接。

【讨论】：

这似乎是我问题的根源。 Excel 如何为多项式拟合与线性回归获得不同的 r 平方值？

你只是给excel一个线性回归的拟合，以及一个多项式模型的拟合吗？它将从两个数据数组中计算 rsq，并假设您从线性模型中对其进行拟合。你给excel什么？什么是 Excel 中的“最适合趋势线”命令？

它是 Excel 绘图功能的一部分。您可以绘制一些数据，右键单击它，然后从几种不同类型的趋势线中进行选择。可以选择查看线的方程以及每种类型的 r 平方值。每种类型的 r 平方值也不同。

@Travis Beale -- 你会为你尝试的每个不同的平均函数得到一个不同的 r 平方（除非两个模型是嵌套的并且较大模型中的额外系数都为 0） .所以当然 Excel 给出了不同的 r 平方值。 @Baltimark - 这是线性回归，所以它是 r 平方的。

【参考方案2】：

根据numpy.polyfit 文档，它拟合线性回归。具体来说，度数为 'd' 的 numpy.polyfit 与均值函数拟合线性回归

E(y|x) = p_d * x**d + p_d-1 * x **(d-1) + ... + p_1 * x + p_0

因此，您只需要计算该拟合的 R 平方即可。 linear regression 上的***页面提供了完整的详细信息。您对 R^2 感兴趣，可以通过多种方式计算，最简单的可能是

SST = Sum(i=1..n) (y_i - y_bar)^2
s-s-reg = Sum(i=1..n) (y_ihat - y_bar)^2
Rsquared = s-s-reg/SST

我使用 'y_bar' 作为 y 的平均值，使用 'y_ihat' 作为每个点的拟合值。

我对 numpy 不是很熟悉（我通常在 R 中工作），所以可能有一种更简洁的方法来计算你的 R 平方，但以下应该是正确的

import numpy

# Polynomial Regression
def polyfit(x, y, degree):
    results = 

    coeffs = numpy.polyfit(x, y, degree)

     # Polynomial Coefficients
    results['polynomial'] = coeffs.tolist()

    # r-squared
    p = numpy.poly1d(coeffs)
    # fit values, and mean
    yhat = p(x)                         # or [p(z) for z in x]
    ybar = numpy.sum(y)/len(y)          # or sum(y)/len(y)
    s-s-reg = numpy.sum((yhat-ybar)**2)   # or sum([ (yihat - ybar)**2 for yihat in yhat])
    sstot = numpy.sum((y - ybar)**2)    # or sum([ (yi - ybar)**2 for yi in y])
    results['determination'] = s-s-reg / sstot

    return results

【讨论】：

我只想指出，使用 numpy 数组函数而不是列表理解会更快，例如numpy.sum((yi - ybar)**2) 更易阅读

根据wiki页面en.wikipedia.org/wiki/Coefficient_of_determination，R^2的最一般定义是R^2 = 1 - SS_err/SS_tot，R^2 = SS_reg/SS_tot只是一个特例。

非线性最小二乘函数的 R 平方怎么样？

【参考方案3】：

r-squareds 上的***文章表明它可以用于一般模型拟合，而不仅仅是线性回归。

【讨论】：

这里很好地描述了 R2 用于非线性回归的问题：blog.minitab.com/blog/adventures-in-statistics/…

【参考方案4】：

一个很晚的回复，但以防万一有人需要为此准备好的功能：

scipy.stats.linregress

即

slope, intercept, r_value, p_value, std_err = scipy.stats.linregress(x, y)

就像@Adam Marples 的回答一样。

【讨论】：

用相关系数分析是合理的，然后做更大的工作，回归。

这个回复只适用于线性回归，是最简单的多项式回归

注意：这里的 r_value 是 Pearson 相关系数，而不是 R 平方。 r_squared = r_value**2

【参考方案5】：

我一直在成功地使用它，其中 x 和 y 是类似数组的。

注意：仅适用于线性回归

def rsquared(x, y):
    """ Return R^2 where x and y are array-like."""

    slope, intercept, r_value, p_value, std_err = scipy.stats.linregress(x, y)
    return r_value**2

【讨论】：

这不是佩拉森的决定系数，而是相关系数的平方——完全是另外一回事。

@liorr 我的理解是决定系数是相关系数的平方

我认为这仅在使用线性回归时才是正确的：en.wikipedia.org/wiki/Coefficient_of_determination“此类情况包括使用 r2 而不是 R2 的简单线性回归。当仅包括截距时，则r2 只是观察到的结果和观察到的预测值之间的样本相关系数（即 r）的平方。”

@liorr 我在我的答案 scipy.stats.linregress 中使用线性回归中的 r**2，所以它是正确的

啊，是的，我没有正确阅读这个问题。在我的辩护中，那是 9 年前的事了，但我还没有。

【参考方案6】：

来自 yanl (yet-another-library) sklearn.metrics 有一个 r2_score 函数；

from sklearn.metrics import r2_score

coefficient_of_dermination = r2_score(y, p(x))

【讨论】：

（注意：“默认值对应于‘variance_weighted’，此行为自0.17版起已弃用，将从0.19开始更改为‘uniform_average’”）

r2_score in sklearn 可能是负值，这不是正常情况。

为什么是r2_score([1,2,3],[4,5,7]) = -16？

我喜欢的一点是它不需要训练模型——我经常从在不同环境中训练的模型计算指标。

【参考方案7】：

我最初发布下面的基准是为了推荐numpy.corrcoef，愚蠢地没有意识到原来的问题已经使用corrcoef，实际上是在询问高阶多项式拟合。我使用 statsmodels 为多项式 r-squared 问题添加了一个实际解决方案，并且我保留了原始基准，这些基准虽然离题，但可能对某人有用。

---

statsmodels 可以直接计算多项式拟合的r^2，这里有两种方法...

import statsmodels.api as sm
import statsmodels.formula.api as smf

# Construct the columns for the different powers of x
def get_r2_statsmodels(x, y, k=1):
    xpoly = np.column_stack([x**i for i in range(k+1)])    
    return sm.OLS(y, xpoly).fit().rsquared

# Use the formula API and construct a formula describing the polynomial
def get_r2_statsmodels_formula(x, y, k=1):
    formula = 'y ~ 1 + ' + ' + '.join('I(x**)'.format(i) for i in range(1, k+1))
    data = 'x': x, 'y': y
    return smf.ols(formula, data).fit().rsquared # or rsquared_adj

要进一步利用statsmodels，还应该查看拟合模型摘要，它可以在 Jupyter/IPython 笔记本中打印或显示为丰富的 HTML 表格。除了 rsquared 之外，结果对象还提供对许多有用统计指标的访问。

model = sm.OLS(y, xpoly)
results = model.fit()
results.summary()

---

以下是我的原始答案，我在其中对各种线性回归 r^2 方法进行了基准测试...

问题中使用的corrcoef 函数仅针对单个线性回归计算相关系数r，因此它不能解决r^2 的问题以进行高阶多项式拟合。但是，不管怎样，我发现对于线性回归，它确实是计算r的最快和最直接的方法。

def get_r2_numpy_corrcoef(x, y):
    return np.corrcoef(x, y)[0, 1]**2

这些是我对 1000 个随机 (x, y) 点的一堆方法进行比较的时间结果：

纯Python（直接r计算） 1000 次循环，3 次中的最佳：每个循环 1.59 毫秒 Numpy polyfit（适用于 n 次多项式拟合） 1000 次循环，3 次中的最佳：每个循环 326 µs Numpy手册（直接r计算） 10000 次循环，3 次中的最佳：每个循环 62.1 µs numpy corrcoef（直接r计算） 10000 次循环，3 次中的最佳：每个循环 56.6 µs Scipy（以r 作为输出的线性回归） 1000 次循环，3 次中的最佳：每个循环 676 µs Statsmodels（可以进行 n 次多项式和许多其他拟合） 1000 次循环，3 次中的最佳：每个循环 422 µs

corrcoef 方法比使用 numpy 方法“手动”计算 r^2 略胜一筹。它比 polyfit 方法快 > 5 倍，比 scipy.linregress 快约 12 倍。只是为了加强 numpy 为您所做的事情，它比纯 python 快 28 倍。我不精通 numba 和 pypy 之类的东西，所以其他人必须填补这些空白，但我认为这对我来说很有说服力，corrcoef 是计算 r 用于简单线性的最佳工具回归。

这是我的基准测试代码。我从 Jupyter Notebook 复制粘贴（很难不称它为 IPython Notebook...），所以如果途中出现任何问题，我深表歉意。 %timeit 魔术命令需要 IPython。

import numpy as np
from scipy import stats
import statsmodels.api as sm
import math

n=1000
x = np.random.rand(1000)*10
x.sort()
y = 10 * x + (5+np.random.randn(1000)*10-5)

x_list = list(x)
y_list = list(y)

def get_r2_numpy(x, y):
    slope, intercept = np.polyfit(x, y, 1)
    r_squared = 1 - (sum((y - (slope * x + intercept))**2) / ((len(y) - 1) * np.var(y, ddof=1)))
    return r_squared
    
def get_r2_scipy(x, y):
    _, _, r_value, _, _ = stats.linregress(x, y)
    return r_value**2
    
def get_r2_statsmodels(x, y):
    return sm.OLS(y, sm.add_constant(x)).fit().rsquared
    
def get_r2_python(x_list, y_list):
    n = len(x_list)
    x_bar = sum(x_list)/n
    y_bar = sum(y_list)/n
    x_std = math.sqrt(sum([(xi-x_bar)**2 for xi in x_list])/(n-1))
    y_std = math.sqrt(sum([(yi-y_bar)**2 for yi in y_list])/(n-1))
    zx = [(xi-x_bar)/x_std for xi in x_list]
    zy = [(yi-y_bar)/y_std for yi in y_list]
    r = sum(zxi*zyi for zxi, zyi in zip(zx, zy))/(n-1)
    return r**2
    
def get_r2_numpy_manual(x, y):
    zx = (x-np.mean(x))/np.std(x, ddof=1)
    zy = (y-np.mean(y))/np.std(y, ddof=1)
    r = np.sum(zx*zy)/(len(x)-1)
    return r**2
    
def get_r2_numpy_corrcoef(x, y):
    return np.corrcoef(x, y)[0, 1]**2
    
print('Python')
%timeit get_r2_python(x_list, y_list)
print('Numpy polyfit')
%timeit get_r2_numpy(x, y)
print('Numpy Manual')
%timeit get_r2_numpy_manual(x, y)
print('Numpy corrcoef')
%timeit get_r2_numpy_corrcoef(x, y)
print('Scipy')
%timeit get_r2_scipy(x, y)
print('Statsmodels')
%timeit get_r2_statsmodels(x, y)

7/28/21 基准测试结果。 （Python 3.7、numpy 1.19、scipy 1.6、statsmodels 0.12）

Python
2.41 ms ± 180 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
Numpy polyfit
318 µs ± 44.3 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
Numpy Manual
79.3 µs ± 4.05 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
Numpy corrcoef
83.8 µs ± 1.37 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
Scipy
221 µs ± 7.12 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
Statsmodels
375 µs ± 3.63 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

【讨论】：

您正在比较 3 种拟合斜率的方法和回归与 3 种不拟合斜率的方法。

是的，我知道这么多......但现在我觉得很傻，因为没有阅读原始问题并看到它已经使用 corrcoef 并且专门针对高阶多项式处理 r^2...... 现在我为不同目的发布我的基准感到愚蠢。哎呀...

我已经使用statsmodels 更新了我对原始问题的解决方案的答案，并为线性回归 r^2 方法的不必要基准测试道歉，我一直认为它很有趣，但离题信息。

我仍然觉得基准测试很有趣，因为我没想到 scipy 的 linregress 会比 statsmodels 更通用的工作。

注意，np.column_stack([x**i for i in range(k+1)]) 可以在 numpy 中使用 x[:,None]**np.arange(k+1) 或使用 numpy 的 vander 函数进行矢量化，这些函数在列中具有相反的顺序。

【参考方案8】：

这是一个用 Python 和 Numpy 计算 加权 r-squared 的函数（大部分代码来自 sklearn）：

from __future__ import division 
import numpy as np

def compute_r2_weighted(y_true, y_pred, weight):
    sse = (weight * (y_true - y_pred) ** 2).sum(axis=0, dtype=np.float64)
    tse = (weight * (y_true - np.average(
        y_true, axis=0, weights=weight)) ** 2).sum(axis=0, dtype=np.float64)
    r2_score = 1 - (sse / tse)
    return r2_score, sse, tse

例子：

from __future__ import print_function, division 
import sklearn.metrics 

def compute_r2_weighted(y_true, y_pred, weight):
    sse = (weight * (y_true - y_pred) ** 2).sum(axis=0, dtype=np.float64)
    tse = (weight * (y_true - np.average(
        y_true, axis=0, weights=weight)) ** 2).sum(axis=0, dtype=np.float64)
    r2_score = 1 - (sse / tse)
    return r2_score, sse, tse    

def compute_r2(y_true, y_predicted):
    sse = sum((y_true - y_predicted)**2)
    tse = (len(y_true) - 1) * np.var(y_true, ddof=1)
    r2_score = 1 - (sse / tse)
    return r2_score, sse, tse

def main():
    '''
    Demonstrate the use of compute_r2_weighted() and checks the results against sklearn
    '''        
    y_true = [3, -0.5, 2, 7]
    y_pred = [2.5, 0.0, 2, 8]
    weight = [1, 5, 1, 2]
    r2_score = sklearn.metrics.r2_score(y_true, y_pred)
    print('r2_score: 0'.format(r2_score))  
    r2_score,_,_ = compute_r2(np.array(y_true), np.array(y_pred))
    print('r2_score: 0'.format(r2_score))
    r2_score = sklearn.metrics.r2_score(y_true, y_pred,weight)
    print('r2_score weighted: 0'.format(r2_score))
    r2_score,_,_ = compute_r2_weighted(np.array(y_true), np.array(y_pred), np.array(weight))
    print('r2_score weighted: 0'.format(r2_score))

if __name__ == "__main__":
    main()
    #cProfile.run('main()') # if you want to do some profiling

输出：

r2_score: 0.9486081370449679
r2_score: 0.9486081370449679
r2_score weighted: 0.9573170731707317
r2_score weighted: 0.9573170731707317

这对应于formula（mirror）：

其中 f_i 是拟合的预测值，y_av 是观测数据的平均值 y_i 是观测数据值。 w_i 是应用于每个数据点的权重，通常 w_i=1。 SSE 是误差平方和，SST 是总平方和。

---

如果有兴趣，R中的代码：https://gist.github.com/dhimmel/588d64a73fa4fef02c8f (mirror)

【参考方案9】：

来自 scipy.stats.linregress 源。他们使用平均平方和方法。

import numpy as np

x = np.array(x)
y = np.array(y)

# average sum of squares:
ssxm, ssxym, ssyxm, ssym = np.cov(x, y, bias=1).flat

r_num = ssxym
r_den = np.sqrt(ssxm * ssym)
r = r_num / r_den

if r_den == 0.0:
    r = 0.0
else:
    r = r_num / r_den

    if r > 1.0:
        r = 1.0
    elif r < -1.0:
        r = -1.0

【参考方案10】：

这是一个非常简单的 python 函数，假设 y 和 y_hat 是熊猫系列，从实际值和预测值计算 R^2：

def r_squared(y, y_hat):
    y_bar = y.mean()
    ss_tot = ((y-y_bar)**2).sum()
    ss_res = ((y-y_hat)**2).sum()
    return 1 - (ss_res/ss_tot)

【讨论】：

这个公式给出的答案与非平凡数据的 numpy 模块不同。这可能是因为 r_squared 是一个优化问题，对最佳拟合线的斜率和偏移有多种解决方案。

上述函数适用于任何模型，线性、非线性、ML 等……它只查看预测值与实际值之间的差异。每个模型通常会创建不同的 R^2。拟合给定模型涉及通过改变模型的参数来最小化 R^2。具有一个自变量和一个因变量的曲线的直线拟合具有唯一解（局部最小值 == 全局最小值）。更复杂的模型，尤其是带有额外自变量的模型，可能有许多局部最小值，而找到全局最小值可能非常困难。

【参考方案11】：

您可以直接执行此代码，这将找到您的多项式，并会找到您的 R 值如果您需要更多解释，可以在下方留言。

from scipy.stats import linregress
import numpy as np

x = np.array([1,2,3,4,5,6])
y = np.array([2,3,5,6,7,8])

p3 = np.polyfit(x,y,3) # 3rd degree polynomial, you can change it to any degree you want
xp = np.linspace(1,6,6)  # 6 means the length of the line
poly_arr = np.polyval(p3,xp)

poly_list = [round(num, 3) for num in list(poly_arr)]
slope, intercept, r_value, p_value, std_err = linregress(x, poly_list)
print(r_value**2)

【参考方案12】：

使用numpy模块（在python3中测试）：

import numpy as np
def linear_regression(x, y): 
    coefs = np.polynomial.polynomial.polyfit(x, y, 1)
    ffit = np.poly1d(coefs)
    m = ffit[0]
    b = ffit[1] 
    eq = 'y = x + '.format(round(m, 3), round(b, 3))
    rsquared = np.corrcoef(x, y)[0, 1]**2
    return rsquared, eq, m, b

rsquared, eq, m, b = linear_regression(x,y)
print(rsquared, m, b)
print(eq)

输出：

0.013378252355751777 0.1316331351105754 0.7928782850418713 
y = 0.132x + 0.793

注意：r² ≠ R² r² 称为“确定系数” R² 是皮尔逊系数的平方

R²，正式合并为 r²，可能是您想要的，因为它是最小二乘拟合，比 r² 的简单分数要好。 Numpy 不怕称其为“corrcoef”，这假定 Pearson 是事实上的相关系数。

【讨论】：

我发布了这个解决方案，因为***文章公式给出的结果与 numpy 解决方案不同。我相信 numpy 模块是正确的，因为 wikipedia 公式不考虑存在多个解决方案（最佳拟合线的不同斜率和偏移量），并且 numpy 显然解决了实际的优化问题，而不仅仅是计算总和的一小部分。 [simple] wikipedia 公式错误的证据是它产生负 r_squared 值，这意味着它为非平凡数据的最佳拟合线提出了错误的斜率。