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求离散序列导数的示例就是小编分享给大家的全部内容了，希望能给大家一个参考，也希望大家多多支持我们。

赞 (0)

相关推荐

  • Python求导数的方法

    本文实例讲述了Python求导数的方法.分享给大家供大家参考.具体实现方法如下: def func(coeff): sum='' for key in coeff: sum=sum+'+'+str(key)+'*'+'x'+'**'+str(coeff[key]) return sum[1:] from sympy import * from sympy.core.sympify import SympifyError expr = func({2:0,3:1,4:2,5:7}) x = Sym

  • python 寻找离散序列极值点的方法

    使用 scipy.signal 的 argrelextrema 函数(API),简单方便 import numpy as np import pylab as pl import matplotlib.pyplot as plt import scipy.signal as signal x=np.array([ 0, 6, 25, 20, 15, 8, 15, 6, 0, 6, 0, -5, -15, -3, 4, 10, 8, 13, 8, 10, 3, 1, 20, 7, 3, 0 ])

  • 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个点拟合成一个函数

  • python 求定积分和不定积分示例

    求f(x) = sin(x)/x 的不定积分和负无穷到正无穷的定积分 sin(x)/x 的不定积分是信号函数sig ,负无穷到正无穷的定积分为pi import math import numpy as np import matplotlib.pyplot as plt from sympy import * #用于求导积分等科学计算 def draw_plot_set():#设置画图格式 ax = plt.gca() #改变坐标轴位置 ax.spines['right'].set_color

  • Python求字符串的长度示例代码

    目录 Python求字符串的长度 补充:Python计算字符串长度的函数 Python三种计算字符串长度的函数 1.使用内置函数len 2.使用for循环 3.使用while循环和切片 Python求字符串的长度 python 写一个函数,求一个字符串的长度,在main函数中输入字符串,并输出其长度. def func(a): #定义一个求字符串长度的函数 '求一个字符串的长度' len=0 for i in a: len+=1 return len if __name__=='__main__

  • 利用python求相邻数的方法示例

    前言 本文主要给大家介绍了关于利用python求相邻数的相关内容,分享出来供大家参考学习,下面话不多说了,来一起看看详细的介绍: 什么是相邻数? 比如5,相邻数为4和6,和5相差1的数,连续相差为1的一组数 需求: 遍历inputList 所有数字,取出所有数字,判断是否有相邻数, 不相邻数字 和 相邻数字 都以 "数组"形式 添加到 outputList 中, 并且 每个"数组" 里 第一位 递减 补全两位数,末位 递增 补全两位数, 每一个数不能小于0, 不能大

  • python开发任意表达式求值全功能示例

    目录 正文 完整的源代码 正文 在之前的基础上进一步实现了全功能表达式求值. 已支持浮点数 已支持字符串的处理,前加一个"(类似lisp语法) 支持减号/负号,一符两用机制 支持所有算术运算符,包括**,//, % 支持全部7个比较运算符 支持与或非3个逻辑运算符 支持自定义数学函数(代码中预设sin函数作为示范) 支持外部提供的变量机制 支持外部设置函数(代码中预设isvar函数作为示范) 支持列表 字典的支持,体现在外部的变量中 结构清晰,易于扩展 具有实用性及学习性 与其说距离DSL只有一

  • Python求一批字符串的最长公共前缀算法示例

    Python求一批字符串的最长公共前缀算法示例

    本文实例讲述了Python求一批字符串的最长公共前缀算法.分享给大家供大家参考,具体如下: 思路一:这个题一拿到手,第一反应就是以第一个字符串strs[0]为标准,如果其他字符串的第一个字符和str[0]的第一个字符串相同,则再比较第二个字符串,以此类推直到出现不同为止. def longestCommonPrefix(self, strs): """ :type strs: List[str] :rtype: str """ if not st

  • 使用python求斐波那契数列中第n个数的值示例代码

    斐波那契数列(Fibonacci sequence),又称黄金分割数列.因数学家列昂纳多·斐波那契(Leonardoda Fibonacci)以兔子繁殖为例子而引入,故又称为"兔子数列",指的是这样一个数列:1.1.2.3.5.8.13.21.34.--在数学上,斐波纳契数列以如下被以递归的方法定义:F(1)=1,F(2)=1, F(n)=F(n-1)+F(n-2)(n>=2,n∈N*) 求斐波那契数列中第n个数的值:1,1,2,3,5,8,13,21,34- 方法一:用for循

  • python opencv之SIFT算法示例

    本文介绍了python opencv之SIFT算法示例,分享给大家,具体如下: 目标: 学习SIFT算法的概念 学习在图像中查找SIFT关键的和描述符 原理: (原理部分自己找了不少文章,内容中有不少自己理解和整理的东西,为了方便快速理解内容和能够快速理解原理,本文尽量不使用数学公式,仅仅使用文字来描述.本文中有很多引用别人文章的内容,仅供个人记录使用,若有错误,请指正出来,万分感谢) 之前的harris算法和Shi-Tomasi 算法,由于算法原理所致,具有旋转不变性,在目标图片发生旋转时依然

  • Python实现梯度下降法的示例代码

    Python实现梯度下降法的示例代码

    目录 1.首先读取数据集 2.初始化相关参数 3.定义计算代价函数–>MSE 4.梯度下降 5.执行 1.首先读取数据集 导包并读取数据,数据自行任意准备,只要有两列,可以分为自变量x和因变量y即可即可. import numpy as np import matplotlib.pyplot as plt data = np.loadtxt("data.csv", delimiter=",") x_data = data[:, 0] y_data = data

随机推荐