Unscented Transformation（アンセンテッド変換）のPython実装

Unscented Transformation（U変換）とは

Unscented Transformation（U変換）は、確率変数 $x$ が非線形関数 $f$ によって変換された確率変数 $y = f(x)$ の平均 $\bar{y}$ と共分散 $P_y$ を推定する手法です。

$x$ の平均 $\bar{x}$ と共分散 $P_x$ が既知である場合、単純にモンテカルロ法で多数のサンプルを生成して変換し、その平均と共分散を計算することも考えられます。

$$ \bar{y} \approx \frac{1}{N} \sum_{i=1}^N f(x_i) $$ $$ P_y \approx \frac{1}{N} \sum_{i=1}^N (f(x_i) - \bar{y})(f(x_i) - \bar{y})^T $$

しかし、この方法では高精度な推定を得るために膨大な数のサンプル $N$ が必要となり、計算コストが非常に高くなります。

U変換は、この問題を解決するために、線形近似を行うことなく、少数の代表点（シグマ点）を用いることで、変換後の確率変数の統計的性質を効率的かつ高精度に推定します。これにより、モンテカルロ法の利点（非線形性への対応）を保ちつつ、計算量を大幅に削減できます。

参考: UKF （Unscented Kalman Filter）って何？

Pythonによる実装

ここでは、2次元の入力 $X=(X_1, X_2)$ を1次元の出力 $Y=f(X)=X_1 \cdot X_2$ に変換する例を考えます。

1. 入力確率変数の設定

入力 $X$ の平均ベクトル $\bar{x}$ と共分散行列 $P_x$ を定義します。

$$ \bar{x} = [E[X_1], E[X_2]] $$ $$ P_x = \begin{pmatrix} \text{Var}[X_1] & \text{Cov}[X_1, X_2] \ \text{Cov}[X_2, X_1] & \text{Var}[X_2] \end{pmatrix} $$

例として、$X_1$ が平均0、分散1、$X_2$ が平均1、分散4、共分散が2であると仮定します。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
import matplotlib.pyplot as plt
import numpy as np
import random
import math
import scipy.linalg

# 入力xの次元数
n = 2
# 出力yの次元数
m = 1

# xの平均ベクトルと共分散行列
x_mean = np.array([0.0, 1.0])
x_P = np.array([[1.0, 2.0], [2.0, 4.0]])

print(f"入力xの平均: {x_mean}")
print(f"入力xの共分散行列:\n{x_P}")

# 非線形変換関数 f(x) = x[0] * x[1]
def f(x):
  return np.array([x[0] * x[1]])

2. シグマ点の計算

シグマ点は、入力確率変数 $x$ の平均と共分散を正確に表現するように選ばれた代表点です。これらの点は、以下の式で計算されます。

$$ \sigma_0 = \bar{x} \tag{1} $$ $$ \sigma_i = \bar{x} + (\sqrt{(n+\lambda)P_x})_i \quad \text{for } i=1, \dots, n \tag{2} $$ $$ \sigma_i = \bar{x} - (\sqrt{(n+\lambda)P_x})_i \quad \text{for } i=n+1, \dots, 2n \tag{3} $$

ここで、$ (\sqrt{(n+\lambda)P_x})_i $ は、行列の平方根（コレスキー分解など）の $i$ 列目を表します。

パラメータ $\lambda$ は、以下の式で計算されます。 $$ \lambda = \alpha^2 (n + \kappa) - n \tag{4} $$

$\alpha$: シグマ点の平均からの広がりを決定するスカラー値（通常0から1の間の小さい正の値）。
$\kappa$: 通常は0に設定される調整パラメータ。シグマ点の広がりをさらに調整します。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
# パラメータ設定
alpha = 0.5
kappa = 0.0

# 式4: lambdaの計算
lam = alpha**2 * (n + kappa) - n
print(f"lambda: {lam}")

# シグマ点を格納する配列 (2n+1個のシグマ点)
sigma_points = np.zeros((n, 2 * n + 1))

# 式1: 最初のシグマ点は平均自身
sigma_points[:, 0] = x_mean

# (n+lambda) * Px の平方根を計算 (コレスキー分解を使用)
sqrt_term = scipy.linalg.cholesky((n + lam) * x_P, lower=True)

# 式2, 3: 残りのシグマ点を計算
for i in range(n):
    sigma_points[:, i + 1] = x_mean + sqrt_term[:, i]
    sigma_points[:, i + n + 1] = x_mean - sqrt_term[:, i]

print("\n生成されたシグマ点:")
for i in range(2 * n + 1):
  print(f"  シグマ点 {i}: {sigma_points[:, i]}")

3. 変換と統計量の計算

計算されたシグマ点それぞれを非線形関数 $f$ で変換し、$y_{\sigma}$ を求めます。 $$ y_{\sigma, i} = f(\sigma_i) \tag{5} $$

次に、各シグマ点に対応する重み $w_i$ を計算します。 $$ w_0 = \frac{\lambda}{n + \lambda} \tag{6} $$ $$ w_i = \frac{1}{2(n + \lambda)} \quad \text{for } i=1, \dots, 2n \tag{7} $$

最後に、これらの重みと変換されたシグマ点 $y_{\sigma}$ を用いて、$y$ の平均 $\bar{y}$ と共分散 $P_y$ を計算します。 $$ \bar{y} \approx \sum_{i=0}^{2n} w_i y_{\sigma, i} \tag{8} $$ $$ P_y \approx \sum_{i=0}^{2n} w_i (y_{\sigma, i} - \bar{y})(y_{\sigma, i} - \bar{y})^T \tag{9} $$

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
# 変換されたシグマ点を格納する配列
sigma_y = np.zeros((m, 2 * n + 1))
# 式5: 非線形変換
for i in range(2 * n + 1):
  sigma_y[:, i] = f(sigma_points[:, i])

# 重み関数を格納する配列
w = np.zeros(2 * n + 1)
# 式6, 7: 重み関数の計算
w[0] = lam / (n + lam)
for i in range(1, 2 * n + 1):
  w[i] = 1 / (2 * (n + lam))

print(f"\n計算された重み: {w}")

# 式8: yの平均の計算
y_mean = np.sum(w * sigma_y, axis=1)

# 式9: yの共分散の計算
y_P = np.zeros((m, m))
for i in range(2 * n + 1):
  diff = sigma_y[:, i] - y_mean
  y_P += w[i] * np.outer(diff, diff) # outer積で(m,m)行列を生成

print(f"\n出力yの推定平均: {y_mean}")
print(f"出力yの推定共分散行列:\n{y_P}")

この実装により、非線形変換後の確率変数の平均と共分散を、少数のシグマ点から効率的に推定できることがわかります。このU変換は、Unscented Kalman Filter (UKF) の基盤となる技術です。