Python求离散序列导数的示例
作者:CaspianR 发布时间:2022-10-04 11:18:52
有一组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阶多项式不能表征的话,更高阶的呢
总体来说,效果还是可以接受的。
如果下阶段找到好的高斯函数拟合方法,会继续更新。
来源:https://blog.csdn.net/renjunsong0/article/details/70176119
猜你喜欢
- 对于SQL的新手,NULL值的概念常常会造成混淆,他们常认为NULL是与空字符串''相同的事。情况并非如此。例如,下述语句是
- 本文实例讲述了php查询mysql数据库并将结果保存到数组的方法。分享给大家供大家参考。具体分析如下:这里主要用到了mysql_fetch_
- 位运算,赋值状态时异或对应位数1的整形,判断状态则与运算对应位数1的整形。最大用处就是同时判断32位状态,节省存储空间,便于扩展, 
- asp fso type属性取得文件类型代码是用来返回类型指定的文件或文件夹。语法FileObject.Type FolderObject.
- 根据我最近的一些实践以及在和一些读者进行关于HTML表格的使用问题沟通之后,决定写这篇文章。总的来说,我注意到由于误导性信息,他们对于tab
- 如何做一个专门显示文本文件的页面? 代码如下:txt.asp<html><head&g
- 前言:看本教程,必须先仔细看前言的内容,否则会进入误区!最近在做个性休闲服装内网站的设计课程,过程中发现,个性元素的应用成为最难的问题,第一
- 前言闲暇时间抽个空写了个三国杀武将手册的小程序,中间有个需求设计的是合成武将皮肤图、竖排的武将姓名、以及小程序码,然后提供保存图片到相册,最
- 有一台windows服务器上跑着mysql的一些应用,现在需要将mysql的数据每天备份,并通过ftp上传到指定的存储服务器上要是在linu
- 天冷,人懒,事多,我就不全文翻译了。只列几个标题,很多内容完全按照我自己的理解写了一下。想读原汁原味的请移步:Icon design tre
- 引言今年互联网的就业环境真的好糟糕啊,好多朋友被优化。我们平常在工作中除了撸好代码,跑通项目之外,还要注意内外兼修。内功和招式都得练👌,才能
- 内容摘要:ASP开发人员为了在他们的设计项目中获得更好的性能和可扩展性而不断努力。幸运地是,有许多书籍和站点在这方面提供了很好的建议。但是这
- 最简单的方法:取整后判断是否和原值相等!javascript的取整函数是:parseIntif(parseInt(value)==value
- 数学函数 1.绝对值 S:select abs(-1) value O:select abs(-1) value from dual 2.取
- 大家应该经常看到在文本框里提示文字,然后一点就没了。通常做法都是默认给个value,通过js来处理。详细实现都不介绍了,大家都会。现在来看一
- <% dim week_ymd(8) '测出可以手动设定日期,比如this_ymd=#2008-04-1
- 如果不清楚字符串的编码格式的话,就可以将这段字符这样检查:$encode = mb_detect_encoding($string, arr
- 问题:1. 访问 ASP 页面时,出现以下错误:Active Server Pages 错误 'ASP 0201'错误无效的
- 索引是快速搜索的关键。MySQL索引的建立对于MySQL的高效运行是很重要的。下面介绍几种常见的MySQL索引类型。在数据库表中,对字段建立
- 今天我去隽辰的博客去看他的文章,在读完他的文章之后,我很自然的就去读网友们给他留的评论,在读的时候我发现他的评论是顺序的,也就是最早的评论在