Python求离散序列导数的示例

yipeiwu_com5年前Python基础

有一组4096长度的数据,需要找到一阶导数从正到负的点,和三阶导数从负到正的点,截取了一小段。

394.0
388.0
389.0
388.0
388.0
392.0
393.0
395.0
395.0
394.0
394.0
390.0
392.0

按照之前所了解的,对离散值求导其实就是求差分,例如第i点的导数(差分)为:

即在一个宽度为2m+1的窗口内通过计算前后m个值加权后的和得到。但是在实际使用过程中效果不是很好。于是想到了同样在一个宽度为2k+1的窗口内,将这2k+1个点拟合成一个函数,然后求导就可以得到任意阶数的导数值。

首先是函数拟合,使用from scipy.optimize import leastsq即最小二乘拟合

from scipy.optimize import leastsq
class search(object):
  def __init__(self, filename):
    self.filename = filename

  def func(self, x, p):
    f = np.poly1d(p)
    return f(x)

  def residuals(self, p, x, y, reg):
    regularization = 0.1 # 正则化系数lambda
    ret = y - self.func(x, p)
    if reg == 1:
      ret = np.append(ret, np.sqrt(regularization) * p)
    return ret

  def LeastSquare(self, data, k=100, order=4, reg=1, show=1): # k为求导窗口宽度,order为多项式阶数,reg为是否正则化
    l = self.len
    step = 2 * k + 1
    p = [1] * order
    for i in range(0, l, step):
      if i + step < l:
        y = data[i:i + step]
        x = np.arange(i, i + step)
      else:
        y = data[i:]
        x = np.arange(i, l)
      try: 
        r = leastsq(self.residuals, p, args=(x, y, reg))
      except:
        print("Error - curve_fit failed")
      fun = np.poly1d(r[0]) # 返回拟合方程系数
      df_1 = np.poly1d.deriv(fun) # 求得导函数
      df_2 = np.poly1d.deriv(df_1)
      df_3 = np.poly1d.deriv(df_2)
      df_value = df_1(x)
      df3_value = df_3(x)

fun = np.poly1d(r[0]),fun返回的是一个 polynomial class,具体使用可以见官方文档numpy.poly1d
polynomial对象可以使用deriv方法求导数,求得的依然是 polynomial对象。 df_value = df_1(x)所得到的就是x这个几个点求得的导数值。

看似大功告成,但是求导的结果并不是很好,如下图,实际最高点在100左右,但是拟合出来的曲线最高点在120左右,而原因在于使用多项式拟合很难准确拟合曲线。

于是想用高斯函数来实现对曲线的拟合,在matlab中试了下,三阶高斯拟合可以很好的拟合曲线,

但是numpy以及sicpy中没有找到类似poly1d这种对象,虽然可以自己定义高斯函数,如下

  def gaussian(self, x, *param):
    fun = param[0]*np.exp(-np.power(x - param[2], 2.) / (2 * np.power(param[4],    2.)))+param[1]*np.exp(-np.power(x - param[3], 2.) / (2 * np.power(param[5], 2.)))
    return fun

但是,在通过最小二乘拟合得到函数参数后只能得到拟合后的点,无法直接求导数..所以并不适合。

所以还是只能回到多项式拟合,如果4阶多项式不能表征的话,更高阶的呢

总体来说,效果还是可以接受的。

如果下阶段找到好的高斯函数拟合方法,会继续更新。

以上这篇Python求离散序列导数的示例就是小编分享给大家的全部内容了,希望能给大家一个参考,也希望大家多多支持【听图阁-专注于Python设计】。

相关文章

Python脚本实现12306火车票查询系统

Python脚本实现12306火车票查询系统

最近我看到看到使用python实现火车票查询,我自己也实现了,感觉收获蛮多的,下面我就把每一步骤都详细给分享出来。(注意使用的是python3) 首先我将最终结果给展示出来: 在cmd命...

Python基于Flask框架配置依赖包信息的项目迁移部署

Python基于Flask框架配置依赖包信息的项目迁移部署

一般在本机上完成基于Flask框架的代码编写后,如果有接口或者数据操作方面需求需要把代码部署到指定服务器上。 一般情况下,使用Flask框架开发者大多数都是选择Python虚拟环境来运行...

Python及Pycharm安装方法图文教程

Python及Pycharm安装方法图文教程

Python及Pycharm安装方法,供大家参考,具体内容如下 1、任务简介 为了学习Python我今天对它进行了安装,整个安装过程忘了截图,故我在虚拟机中重新安装了一遍,并将Pytho...

python3.8与pyinstaller冲突问题的快速解决方法

python3.8与pyinstaller冲突问题的快速解决方法

安装pyinstaller 安装的时候 进入cmd pip install pyinstaller 发现安装报错! 解决办法: # 自主下载pyinstaller包,进行手动安装...

Python的collections模块中namedtuple结构使用示例

namedtuple 就是命名的 tuple,比较像 C 语言中 struct。一般情况下的 tuple 是 (item1, item2, item3,...),所有的 item 都只能...