カルマンフィルタとベイズ推定によるノイズ付きサイン波


## 🎯 はじめに

現実世界のデータ（音、加速度、信号など）は常に**ノイズ**を含みます。
このノイズを抑えて「真の信号」を推定するには、\*\*カルマンフィルタ（Kalman Filter）**や**ベイズ推定（Bayesian Inference）\*\*といった手法が非常に有効です。

本稿では、以下の手法を統合した高度な信号復元方法を紹介します：

* サイン波（正弦波）にノイズを加えた時系列データを生成
* **調和振動子（harmonic oscillator）モデル**を状態遷移行列として利用
* **カルマンフィルタによる動的推定**
* **ベイズ推定による確率的な信念更新**
* 両者を可視化して比較

---

## 📐 背景理論

### ◼️ ノイズ付きサイン波の生成

真の信号：

$$
x_t = \sin(2\pi f t)
$$

観測値（ノイズ付き）：

$$
z_t = x_t + \epsilon_t,\quad \epsilon_t \sim \mathcal{N}(0, \sigma^2)
$$

---

## 🤖 モデル構築：状態空間表現

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4046871/a373efe8-ba9c-4838-9cec-2120624ad40c.png)




---

### ◼️ 観測行列（H）

我々の観測は状態の「位置成分」だけです：

$$
H = \begin{bmatrix} 1 & 0 \end{bmatrix}
$$

---
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4046871/7f81dd50-fcfe-4d74-8e83-46f2d5c10e79.png)

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4046871/1544354f-0b5d-4672-ae44-8fe23855cb68.png)

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4046871/842959e6-9664-41ad-a738-b270f163187a.png)



---

###  ベイズ推定（Bayesian Update）

カルマン更新のうち「位置成分」だけを切り出して、ベイズの定理を単独で適用する：

* 事前分布：

  $$
  p(x) = \mathcal{N}(\mu_{\text{prior}}, \sigma^2_{\text{prior}})
  $$
* 尤度：

  $$
  p(z|x) = \mathcal{N}(z, \sigma^2_{\text{obs}})
  $$
* 事後分布：

  $$
  p(x|z) = \mathcal{N}(\mu_{\text{post}}, \sigma^2_{\text{post}})
  $$

#### ◼️ 事後分布のパラメータ：

* 分散：

  $$
  \sigma^2_{\text{post}} = \left(\frac{1}{\sigma^2_{\text{prior}}} + \frac{1}{\sigma^2_{\text{obs}}}\right)^{-1}
  $$
* 平均：

  $$
  \mu_{\text{post}} = \sigma^2_{\text{post}} \left(\frac{\mu_{\text{prior}}}{\sigma^2_{\text{prior}}} + \frac{z}{\sigma^2_{\text{obs}}} \right)
  $$

---

## 📊 結果と可視化の解釈

### 📈 表示されるプロット

| 色   | 意味                |
| --- | ----------------- |
| 青線  | 真のサイン波            |
| 灰色点 | ノイズ付き観測値          |
| 橙線  | カルマンフィルタによる推定（動的） |
| 黒破線 | ベイズ推定による平均（静的）    |
| 灰帯  | ベイズ推定の1σ信頼区間      |

### 🔍 観察ポイント

* カルマン推定は時間的文脈を活かすためノイズに強い
* ベイズ推定は1点ごとの確率的信念を示し、**信頼区間の可視化**に優れる
* 両者の併用により、**動的かつ不確かさを伴う信号復元**が実現

---



```python
# Program Name: kalman_sine_bayes.py
# Creation Date: 20250603
# Overview: High-accuracy Kalman Filter + Bayesian update for noisy sine wave estimation
# Usage: Run in Python environment with required libraries installed

# -------------------- Install Required Libraries --------------------
!pip install numpy matplotlib scipy

# -------------------- Import Libraries --------------------
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm

# -------------------- Parameters (一元管理) --------------------
np.random.seed(42)             # 乱数シード / Random seed
n_steps = 100                  # 時系列ステップ数 / Number of time steps
true_freq = 0.2                # サイン波の周波数 / Frequency of sine wave
omega = 2 * np.pi * true_freq # 角周波数 / Angular frequency
noise_std = 0.5                # 観測ノイズ標準偏差 / Observation noise std
process_var = 1e-6             # システムノイズ分散 / Process noise variance
obs_var = noise_std**2         # 観測ノイズ分散 / Observation noise variance
dt = 1.0                       # 時間間隔 / Time interval

# -------------------- Generate Noisy Sine Wave / ノイズ付きサイン波の生成 --------------------
t = np.arange(n_steps)
true_signal = np.sin(omega * t)                                        # 真のサイン波 / True sine wave
observations = true_signal + np.random.normal(0, noise_std, n_steps)  # ノイズ付き観測値 / Noisy observations

# -------------------- Kalman Filter Initialization / カルマンフィルタの初期化 --------------------
x_est = np.array([0.0, omega])      # 初期状態 [信号, 速度] / Initial state: [position, velocity]
P = np.eye(2) * 1.0                 # 共分散行列 / Initial covariance matrix

# 状態遷移行列 (調和振動子モデル) / State transition matrix (harmonic oscillator)
A = np.array([
    [np.cos(omega * dt), np.sin(omega * dt) / omega],
    [-omega * np.sin(omega * dt), np.cos(omega * dt)]
])

H = np.array([[1, 0]])              # 観測行列 / Observation model
Q = process_var * np.eye(2)         # システムノイズ共分散 / Process noise covariance
R = np.array([[obs_var]])           # 観測ノイズ共分散 / Observation noise covariance

# -------------------- Storage Arrays / 推定結果記録用リスト --------------------
kalman_estimates = []              # カルマン推定値 / Kalman estimates
bayes_means = []                   # ベイズ推定平均 / Bayesian posterior means
bayes_vars = []                    # ベイズ推定分散 / Bayesian posterior variances

# -------------------- Kalman Filter + Bayesian Estimation Loop --------------------
for z in observations:
    # ---- Prediction Step / 予測ステップ ----
    x_pred = A @ x_est
    P_pred = A @ P @ A.T + Q

    # ---- Kalman Update Step / カルマン更新ステップ ----
    y = z - float(H @ x_pred)                              # 残差 / Innovation
    S = float(H @ P_pred @ H.T + R)                        # 残差の分散 / Innovation covariance
    K = P_pred @ H.T / S                                   # カルマンゲイン / Kalman gain
    x_est = x_pred + K.flatten() * y                       # 状態更新 / Updated state estimate
    P = (np.eye(2) - K @ H) @ P_pred                       # 分散更新 / Updated covariance
    kalman_estimates.append(float(x_est[0]))               # 位置成分記録 / Store position

    # ---- Bayesian Update Step / ベイズ更新ステップ ----
    prior_mean = float(x_pred[0])
    prior_var = float(P_pred[0, 0])
    post_var = 1 / (1 / prior_var + 1 / obs_var)           # 事後分散 / Posterior variance
    post_mean = post_var * (prior_mean / prior_var + z / obs_var)  # 事後平均 / Posterior mean

    bayes_means.append(post_mean)
    bayes_vars.append(post_var)

# -------------------- Plot Results / 結果の描画 --------------------
plt.figure(figsize=(12, 6))
plt.plot(t, true_signal, label='True Sine Wave')
plt.plot(t, observations, label='Noisy Observations', alpha=0.4)
plt.plot(t, kalman_estimates, label='Kalman Filter Estimates')
plt.plot(t, bayes_means, label='Bayesian Posterior Means', linestyle='--')
plt.fill_between(t,
                 np.array(bayes_means) - np.sqrt(bayes_vars),
                 np.array(bayes_means) + np.sqrt(bayes_vars),
                 color='gray', alpha=0.3, label='Bayesian 1σ Interval')
plt.xlabel('Time step')
plt.ylabel('Signal')
plt.title('Kalman + Bayesian Estimation of Noisy Sine Wave')
plt.legend()
plt.grid(True)
plt.show()


```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4046871/99f67f23-8a64-4916-b496-bd68e30a63ae.png)