请选择 进入手机版 | 继续访问电脑版
爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 18855|回复: 11

[源代码] DIY滑动平均算法,解决了端点问题,Python版

[复制链接]

新浪微博达人勋

发表于 2021-9-27 17:26:43 | 显示全部楼层 |阅读模式

登录后查看更多精彩内容~

您需要 登录 才可以下载或查看,没有帐号?立即注册 新浪微博登陆

x
本帖最后由 15195775117 于 2021-9-27 17:28 编辑

今天处理数据才突然发现,我一直没有掌握一种良好的滑动平均算法。

滑动平均的计算原理,导致其存在端点问题,

以前用IDL语言的smooth函数,也常感觉两端的数怪怪的,并未深究,


今天研究了一下,发现smooth函数其实只是将端点按原值连到均值序列上而已!

1310239080.jpg

网上还到处流传一个Python版的滑动平均处理方法:
https://blog.csdn.net/weixin_42232219/article/details/90328880


我将其提供的3种处理模式一一测试,发现该算法的端点处理效果十分不好:
292809199.jpg 2130384586.jpg 211720656.jpg


所以我自己写了一个算法:


'''
<该函数实现滑动平均>
滑动平均是这样:
一个100个元素的序列,从第一个数开始,每5个数字做一次平均,
但开头2个数和最后2个数的位置将会没有对应的均值,真正计算出来的均值是中间的96个数,
这就存在了端点问题。

我打算这样处理端点的值:
既然端点与边界的距离,已经小于步长的半径,那就将它与边界的距离算它的步长,
例如步长=9,那半径=4,从头数位置第3的数,虽然不能按步长9来取平均,
但可以按步长5来取平均,从头数位置第2的数,可以按步长3来取平均.
对于两端的2个数,将头尾两个数的均值算成是两个端点位置的均值,这样看起来比较自然.

扩展:
1、求均值可换成求中值
2、1维平滑可扩展为2维平滑

'''
import numpy as np

# 注意:x是一维序列,r是一个数字

def smooth(x,r):
    n=len(x)
    #准备一个空序列用来装中间可计算均值的部分的计算结果:
    body=[]

    # 从第r到第n-r的位置(不包括最后一个你懂的),
    # 将每个位置前后总计(2r+1)个数求平均,追加到body中:
    for i in range(r,n-r):
        body.append(np.mean(x[i-r:i+r+1]))

    #准备一个空序列用来装头部边沿的计算结果:
    head=[]
    for i in range(0,r):
        # 该位置两端最大的空间半径:
        r2=i
        # 如果该位置不是端点:
        if r2>0:
            zone=x[i-r2:i+r2+1]
            head.append(np.mean(zone))
        # 如果该位置是端点:
        if r2==0:
            head.append(np.mean(x[0:2]))

    #准备一个空序列用来装尾部边沿的计算结果:
    tail=[]
    for i in range(n-r,n):
        # 该位置两端最大的空间半径:
        r2=n-i-1
        # 如果该位置不是端点:
        if r2>0:
            zone=x[i-r2:i+r2+1]
            tail.append(np.mean(zone))
        # 如果该位置是端点:
        if r2==0:
            tail.append(np.mean(x[-2:]))

    # 把头,身,尾连起来:
    head.extend(body)
    head.extend(tail)
    return head

if __name__=='__main__':
    import matplotlib.pyplot as plt
    # 生成一个n个正态分布随机数的序列,元素值都在[0,100]之间
    n=50
    x=np.random.normal(10,1,(n))
    # 现在取一个滑动平均的步长,
    # (一般这种计算,例如还有滑动取中值,步长最好是奇数,
    # 因为奇数中间的位置好跟原数组对应)
    # 这里我们使用步长的半径,例如步长9,那半径就是4
    r=4
    y=smooth(x,r)
    plt.plot(x,c='blue',marker='o')
    plt.plot(y,c='red',marker='*')
    plt.show()



效果图:
3.png 2.png 1.png



292809199.jpg
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2021-10-2 13:19:43 | 显示全部楼层
厉害了,感谢分享!
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2021-10-2 22:30:35 | 显示全部楼层
river 发表于 2021-10-2 13:19
厉害了,感谢分享!

分享心得,共同进步!
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2021-10-3 09:01:01 | 显示全部楼层
感谢分享,给楼主点赞
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2021-10-4 12:43:59 | 显示全部楼层
scipy.ndimage.gaussian_filter1d、scipy.ndimage.gaussian_filter不香么?
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2021-10-8 08:48:59 | 显示全部楼层
Masterpiece 发表于 2021-10-4 12:43
scipy.ndimage.gaussian_filter1d、scipy.ndimage.gaussian_filter不香么?

灰常感谢!
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2021-10-10 15:27:56 | 显示全部楼层
感谢5楼朋友的良策,这个一维高斯平滑挺好使的!

也不需要考虑边界问题了,返回值与输入值大小相同。

待我有空再试试
gaussian_filter2d


from scipy.ndimage import gaussian_filter1d
import numpy as np
import matplotlib.pyplot as plt

rng = np.random.default_rng()
x = rng.standard_normal(101).cumsum()

y3 = gaussian_filter1d(x, 3)
y6 = gaussian_filter1d(x, 6)

plt.plot(x, 'k', label='original data')
plt.plot(y3, '--', label='filtered, sigma=3')
plt.plot(y6, ':', label='filtered, sigma=6')
plt.legend()
plt.grid()
plt.show()


Figure_1.png


密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2021-10-11 11:05:29 | 显示全部楼层
Masterpiece 发表于 2021-10-4 12:43
scipy.ndimage.gaussian_filter1d、scipy.ndimage.gaussian_filter不香么?

有用,感谢分享
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2022-7-9 09:59:34 | 显示全部楼层
请问一下,用楼主的滑动平均程序,如果x是一维序列,序列当中有缺测值,还会影响结果吗?应该怎么办呢?
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2022-7-9 12:27:32 | 显示全部楼层
qq469015280 发表于 2022-7-9 09:59
请问一下,用楼主的滑动平均程序,如果x是一维序列,序列当中有缺测值,还会影响结果吗?应该怎么办呢?

如有有缺测值,就赋np.nan,然后里面的np.mean()换成np.nanmean(),就行了
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表