快捷導(dǎo)航

Python求離散序列導(dǎo)數(shù)的示例

更新時(shí)間：2019年07月10日 09:33:25 作者：CaspianR

今天小編就為大家分享一篇Python求離散序列導(dǎo)數(shù)的示例，具有很好的參考價(jià)值，希望對(duì)大家有所幫助。一起跟隨小編過(guò)來(lái)看看吧

有一組4096長(zhǎng)度的數(shù)據(jù)，需要找到一階導(dǎo)數(shù)從正到負(fù)的點(diǎn)，和三階導(dǎo)數(shù)從負(fù)到正的點(diǎn)，截取了一小段。

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

按照之前所了解的，對(duì)離散值求導(dǎo)其實(shí)就是求差分，例如第i點(diǎn)的導(dǎo)數(shù)（差分）為：

即在一個(gè)寬度為2m+1的窗口內(nèi)通過(guò)計(jì)算前后m個(gè)值加權(quán)后的和得到。但是在實(shí)際使用過(guò)程中效果不是很好。于是想到了同樣在一個(gè)寬度為2k+1的窗口內(nèi)，將這2k+1個(gè)點(diǎn)擬合成一個(gè)函數(shù)，然后求導(dǎo)就可以得到任意階數(shù)的導(dǎo)數(shù)值。

首先是函數(shù)擬合，使用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 # 正則化系數(shù)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為求導(dǎo)窗口寬度,order為多項(xiàng)式階數(shù),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]) # 返回?cái)M合方程系數(shù)
      df_1 = np.poly1d.deriv(fun) # 求得導(dǎo)函數(shù)
      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返回的是一個(gè) polynomial class，具體使用可以見(jiàn)官方文檔numpy.poly1d
polynomial對(duì)象可以使用deriv方法求導(dǎo)數(shù)，求得的依然是 polynomial對(duì)象。 df_value = df_1(x)所得到的就是x這個(gè)幾個(gè)點(diǎn)求得的導(dǎo)數(shù)值。

看似大功告成，但是求導(dǎo)的結(jié)果并不是很好，如下圖，實(shí)際最高點(diǎn)在100左右，但是擬合出來(lái)的曲線最高點(diǎn)在120左右，而原因在于使用多項(xiàng)式擬合很難準(zhǔn)確擬合曲線。

于是想用高斯函數(shù)來(lái)實(shí)現(xiàn)對(duì)曲線的擬合，在matlab中試了下，三階高斯擬合可以很好的擬合曲線，

但是numpy以及sicpy中沒(méi)有找到類似poly1d這種對(duì)象，雖然可以自己定義高斯函數(shù)，如下

  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

但是，在通過(guò)最小二乘擬合得到函數(shù)參數(shù)后只能得到擬合后的點(diǎn)，無(wú)法直接求導(dǎo)數(shù)..所以并不適合。

所以還是只能回到多項(xiàng)式擬合，如果4階多項(xiàng)式不能表征的話，更高階的呢