kalman滤波python实现——基于维纳退化模型(1)

（1）Python所有方向的学习路线（新版）

这是我花了几天的时间去把Python所有方向的技术点做的整理，形成各个领域的知识点汇总，它的用处就在于，你可以按照上面的知识点去找对应的学习资源，保证自己学得较为全面。

最近我才对这些路线做了一下新的更新，知识体系更全面了。

（2）Python学习视频

包含了Python入门、爬虫、数据分析和web开发的学习视频，总共100多个，虽然没有那么全面，但是对于入门来说是没问题的，学完这些之后，你可以按照我上面的学习路线去网上找其他的知识资源进行进阶。

（3）100多个练手项目

我们在看视频学习的时候，不能光动眼动脑不动手，比较科学的学习方法是在理解之后运用它们，这时候练手项目就很适合了，只是里面的项目比较多，水平也是参差不齐，大家可以挑自己能做的项目去练练。

网上学习资料一大堆，但如果学到的知识不成体系，遇到问题时只是浅尝辄止，不再深入研究，那么很难做到真正的技术提升。

需要这份系统化学习资料的朋友，可以戳这里获取

一个人可以走的很快，但一群人才能走的更远！不论你是正从事IT行业的老鸟或是对IT行业感兴趣的新人，都欢迎加入我们的的圈子（技术交流、学习资源、职场吐槽、大厂内推、面试辅导），让我们一起学习成长！

• 代码实现 ：

import matplotlib.pyplot as plt
from math import pi
import numpy as np
from numpy import log, exp, sqrt
from numpy.random import normal, randn

def kalman_filter_for_weiner_process(t, X):
    ''''X is an array'''
    n_samples = len(X)
    _size = n_samples + 5
    sigma2 = 10                           # difussion coef
    Q = 10                                # drift noise variation
    K = np.zeros(_size)
    P = np.ones([_size,_size])  # P0 = 1

    λ_hat = np.zeros(_size) # a0 = 1
    y_pred = np.zeros(_size)

    for i in np.arange(0, n_samples):
        if i == 0:
            y_pred[i] = X[i]
        else:    
           # Prediction (expectation)
            P[i][i-1] = P[i-1][i-1] + Q 
            K[i] = (t[i] - t[i-1])**2 * P[i][i-1] + sigma2 * (t[i] - t[i-1]) # kalman 增益
            λ_hat[i] = λ_hat[i-1] + P[i][i-1] * (t[i] - t[i-1]) * (X[i] - X[i-1] - λ_hat[i-1] * (t[i] - t[i-1])) / K[i]

            # Correction (variance)
            P[i][i] = P[i][i-1] - P[i][i-1] * (t[i] - t[i-1])**2 / K[i] * P[i][i-1]

            # 基于更新后的λ_hat进行下一步的预测
            y_pred[i] = y_pred[i-1] + λ_hat[i] * (t[i] - t[i-1])
            
        print(f"λ_hat：{λ_hat[i]}，P：{P[i][i]}")

    plt.scatter(t,X, color = 'blue')
    plt.plot(t,X, color = 'blue', label = 'actual')
    plt.plot(t, y_pred[:len(t)], color = 'red', label='kalman')
    plt.scatter(t, y_pred[:len(t)], color = 'red')
    print(y_pred)
    print(λ_hat)
    print(X)
    plt.legend()
    

# data
t = np.arange(1,15)
X = np.sin(np.arange(1,15))  * 100
kalman_filter_for_weiner_process(t, X)

• 输出结果：​​​​​​​​​​​​​​

改进的kalman滤波

• 算法细节：

• 代码实现：

import matplotlib.pyplot as plt
from math import pi
import numpy as np
from numpy import log, exp, sqrt
from numpy.random import normal, randn

def strong_tracking_kalman_filter_for_weiner_process(t, X):
    ''''X is an array'''
    n_samples = len(t)
    _size = n_samples + 5 
    print("n_samples:", n_samples)
    sigma2 = 10                           # difussion coef
    Q = 10                               # drift noise variation
    alpha, rho = 0.5, 0.5
    V = np.zeros(_size)
    gamma = np.zeros(_size)
    B = np.zeros(_size)
    C = np.zeros(_size)
    v = np.zeros(_size)
    

    K = np.zeros(_size)
    P = np.ones([_size,_size])  # P0 = 1

    λ_hat = np.zeros(_size) # a0 = 0
    y_pred = np.zeros(_size)

    for i in np.arange(0, n_samples):
        if i == 0:
            y_pred[i] = X[i]
        else:    
            gamma[i] = X[i] - X[i-1] - λ_hat[i-1] * (t[i] - t[i-1])
            if i == 1:
                V[i] = gamma[i]**2
            elif i > 1:
                V[i] = (rho * V[i-1] + gamma[i]**2) / (1 + rho)
            B[i] = V[i] - Q * (t[i] - t[i-1])**2 - alpha * sigma2 * (t[i] - t[i-1])
            C[i] = P[i-1][i-1] * (t[i] - t[i-1])**2 - alpha * sigma2 * (t[i] - t[i-1])
            v0 = B[i] / C[i]
            if v0 >= 1:
                v[i] = v0
            else:
                v[i] = 1

            # Prediction (expectation)
            P[i][i-1] = v[i] * P[i-1][i-1] + Q 
            K[i] = (t[i] - t[i-1])**2 * P[i][i-1] + sigma2 * (t[i] - t[i-1]) # kalman 增益
            λ_hat[i] = λ_hat[i-1] + P[i][i-1] * (t[i] - t[i-1]) / K[i] * (X[i] - X[i-1] - λ_hat[i-1] * (t[i] - t[i-1])) 

            # Correction (variance)
            P[i][i] = P[i][i-1] - P[i][i-1] * (t[i] - t[i-1])**2 / K[i] * P[i][i-1]

            # 基于调整过后的λ_hat进行下一步的预测
            y_pred[i] = y_pred[i-1] + λ_hat[i] * (t[i] - t[i-1])

        print(f"λ_hat：{λ_hat[i-1]}，P：{P[i][i]}")

    plt.plot(t, X, color = 'blue', label = 'actual')
    plt.scatter(t, X, color = 'blue')


### 最后

> **🍅 硬核资料**：关注即可领取PPT模板、简历模板、行业经典书籍PDF。  
> **🍅 技术互助**：技术群大佬指点迷津，你的问题可能不是问题，求资源在群里喊一声。  
> **🍅 面试题库**：由技术群里的小伙伴们共同投稿，热乎的大厂面试真题，持续更新中。  
> **🍅 知识体系**：含编程语言、算法、大数据生态圈组件（Mysql、Hive、Spark、Flink）、数据仓库、Python、前端等等。




**网上学习资料一大堆，但如果学到的知识不成体系，遇到问题时只是浅尝辄止，不再深入研究，那么很难做到真正的技术提升。**

**[需要这份系统化学习资料的朋友，可以戳这里获取](https://bbs.csdn.net/topics/618317507)**

**一个人可以走的很快，但一群人才能走的更远！不论你是正从事IT行业的老鸟或是对IT行业感兴趣的新人，都欢迎加入我们的的圈子（技术交流、学习资源、职场吐槽、大厂内推、面试辅导），让我们一起学习成长！**