Unscented Transformation（U変換）の理論とPython実装：シグマ点生成・非線形伝播・平均共分散推定（numpy.linalg.cholesky）

numpy.linalg.cholesky と numpy.dot で Unscented 変換を Python 実装。シグマ点生成（2n+1点）・重み係数 Wm/Wc・非線形関数による伝播・変換後の平均と共分散の推定アルゴリズム、なぜ2次モーメントまで厳密に再現されるかの証明、モンテカルロ法との数値比較、UKFへの応用、EKFとの比較、van der Merwe スケーリングまで解説。

19 min read

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であると仮定します。

import numpy as np
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]])

実行結果:

入力xの平均: [0. 1.]
入力xの共分散行列:
[[1. 2.]
 [2. 4.]]

なお、この共分散行列の相関係数を計算すると \(\rho = \mathrm{Cov}(X_1,X_2)/\sqrt{\mathrm{Var}[X_1]\mathrm{Var}[X_2]} = 2/(1\cdot 2) = 1\) となり、\(X_1\) と \(X_2\) が完全に線形従属な縮退したケースになっています。後述するように、これは行列の平方根計算において注意が必要な特異ケースです。

2. シグマ点の計算

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

\[ \sigma_0 = \bar{x} \tag{1} \] \[ \sigma_i = \bar{x} + \left(\sqrt{(n+\lambda) P_x}\right)_i, \quad i=1, \dots, n \tag{2} \] \[ \sigma_i = \bar{x} - \left(\sqrt{(n+\lambda) P_x}\right)_i, \quad i=n+1, \dots, 2n \tag{3} \]

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

パラメータ \(\lambda\) は、以下の式で計算されます。

\[ \lambda = \alpha^2 (n + \kappa) - n \tag{4} \]
  • \(\alpha\) : シグマ点の平均からの広がりを決定するスカラー値（通常0から1の間の小さい正の値）。
  • \(\kappa\) : 通常は0に設定される調整パラメータ。シグマ点の広がりをさらに調整します。

なぜこの重みで平均と共分散が厳密に再現されるのか（証明）

式(1)〜(3)のシグマ点と後述の式(6)(7)の重みを天下り的に受け入れるのではなく、なぜこの \(2n+1\) 個の点と重みの組が、どんな確率分布であっても（ガウス分布を仮定せずとも）入力の平均 \(\bar x\) と共分散 \(P_x\) を厳密に再現するのかを証明します。

\(S=\sqrt{(n+\lambda)P_x}\) を行列平方根（\(SS^T=(n+\lambda)P_x\) ）とし、その第 \(i\) 列を \(s_i\) と書くと、式(2)(3)は \(\sigma_i=\bar x + s_i\) （\(i=1,\dots,n\) ）、\(\sigma_{i+n}=\bar x - s_i\) （\(i=1,\dots,n\) ）と書き直せます。

(a) 重みの総和が1になること

\[ \sum_{i=0}^{2n} w_i = \frac{\lambda}{n+\lambda} + 2n \cdot \frac{1}{2(n+\lambda)} = \frac{\lambda+n}{n+\lambda} = 1 \]

(b) 平均が厳密に再現されること

\[ \sum_{i=0}^{2n} w_i \sigma_i = w_0\bar x + \sum_{i=1}^n \frac{1}{2(n+\lambda)}(\bar x+s_i) + \sum_{i=1}^n \frac{1}{2(n+\lambda)}(\bar x-s_i) \]

\(+s_i\) と \(-s_i\) の項が打ち消し合うため、

\[ \sum_{i=0}^{2n} w_i \sigma_i = \left(w_0+\frac{n}{n+\lambda}\right)\bar x = \left(\frac{\lambda}{n+\lambda}+\frac{n}{n+\lambda}\right)\bar x = \bar x \]

シグマ点が平均を中心に符号対称（\(+s_i\) と \(-s_i\) ）に配置されているため、\(\lambda\) の値や \(P_x\) の中身によらず常に平均が厳密に一致します。

(c) 共分散が厳密に再現されること

\(\sigma_0-\bar x=0\) なので \(w_0\) の項は寄与せず、

\[ \sum_{i=0}^{2n} w_i(\sigma_i-\bar x)(\sigma_i-\bar x)^T = \frac{1}{2(n+\lambda)}\sum_{i=1}^n \left[s_is_i^T + (-s_i)(-s_i)^T\right] = \frac{1}{n+\lambda}\sum_{i=1}^n s_is_i^T \]

行列積の定義より \(\sum_{i=1}^n s_is_i^T = SS^T\) なので、

\[ \sum_{i=0}^{2n} w_i(\sigma_i-\bar x)(\sigma_i-\bar x)^T = \frac{SS^T}{n+\lambda} = \frac{(n+\lambda)P_x}{n+\lambda} = P_x \]

以上(a)〜(c)により、\(n+\lambda \neq 0\) である限り、シグマ点集合はどんな確率分布であっても、その平均・共分散という最初の2つのモーメントを厳密に再現することが証明できました。ガウス分布を一切仮定していない点が重要で、シグマ点は「分布を近似するサンプル点群」ではなく、「指定された平均・共分散を持つ、たった \(2n+1\) 個のアトムからなる離散分布」そのものです（ 信号処理におけるフィルタリング手法の基礎 ではスカラー変数を例にテイラー展開による直感的な説明にとどめていますが、本記事では一般の次元 \(n\) について厳密に証明しています）。

非線形変換後の平均が2次まで正確な理由

出力 \(y=f(x)\) （ここではスカラー出力とします）についても、なぜ2次のテイラー項まで正確に再現できるのかを見ておきます。\(f\) を \(\bar x\) の周りでテイラー展開すると、

\[ f(x) = f(\bar x) + \nabla f(\bar x)^T(x-\bar x) + \frac{1}{2}(x-\bar x)^T H(x-\bar x) + R_3(x) \]

（\(H\) は \(f\) のヘッシアン、\(R_3\) は3次以降の残余項）。シグマ点による加重平均をとると、(b)より1次項の係数 \(\sum_i w_i(\sigma_i-\bar x)=0\) 、(c)より2次項は \(\frac{1}{2}\mathrm{tr}(HP_x)\) に厳密に一致するため、

\[ \sum_{i=0}^{2n} w_i f(\sigma_i) = f(\bar x) + \frac{1}{2}\mathrm{tr}(HP_x) + \sum_{i=0}^{2n} w_i R_3(\sigma_i) \]

一方、真の期待値も同じテイラー展開を確率変数 \(X\) について期待値をとるだけなので、

\[ E[f(X)] = f(\bar x) + \frac{1}{2}\mathrm{tr}(HP_x) + E[R_3(X)] \]

両者の差は3次以降の残余項だけに現れます。特に \(f\) が2次多項式（本記事の \(f(x)=x_1x_2\) など）であれば \(R_3\equiv 0\) が恒等的に成り立つため、シグマ点による平均推定は入力分布の形状によらず厳密に真の平均と一致します。一方、ヤコビ行列だけで線形化するEKFは1次の項までしか捉えないため、この \(\frac{1}{2}\mathrm{tr}(HP_x)\) という2次のバイアス項を完全に見落とします（EKFの場合、伝播後の平均推定値は単に \(f(\bar x)\) です）。後ほど数値例でこの差を具体的に確認します。

実装：行列の平方根とその落とし穴（特異な共分散行列）

式(2)(3)の平方根 \(\sqrt{(n+\lambda)P_x}\) は通常コレスキー分解で計算しますが、本記事の \(P_x=\begin{pmatrix}1&2\\2&4\end{pmatrix}\) は前述の通り相関係数1の特異行列（\(\det P_x=0\) 、固有値は \(\{0,5\}\) ）です。正定値性を要求するコレスキー分解は、このような行列に対して失敗します。

try:
    scipy.linalg.cholesky(x_P, lower=True)
except Exception as e:
    print(f"cholesky failed: {type(e).__name__}: {e}")

実行結果:

cholesky failed: LinAlgError: Internal potrf return info = [2] for slices [0].

このように、共分散行列が特異（あるいは数値誤差でわずかに半正定値を割り込む）場合、コレスキー分解ベースの実装はそのまま例外で停止してしまいます。完全相関という極端なケースに限らず、次元縮退したプロセスノイズや長時間の伝播後の丸め誤差でも同様の問題が起こり得るため、実務では正定値性を要求しないロバストな平方根計算に置き換える必要があります。固有値分解を使えば、半正定値行列（固有値が0以上）に対しても常に平方根を計算できます。

def sym_sqrt(P):
    """半正定値対称行列の平方根を固有値分解で計算する（特異行列でも失敗しない）"""
    eigvals, eigvecs = np.linalg.eigh(P)
    eigvals = np.clip(eigvals, 0, None)  # 数値誤差による微小な負値を0に丸める
    return eigvecs @ np.diag(np.sqrt(eigvals)) @ eigvecs.T


S_check = sym_sqrt(x_P)
print(f"S_check @ S_check.T:\n{S_check @ S_check.T}")

実行結果:

S_check @ S_check.T:
[[1. 2.]
 [2. 4.]]

\(SS^T\) が元の \(P_x\) を機械精度で正確に再現していることが確認できます。以降のシグマ点計算では、この sym_sqrt を使用します。

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} \]

パラメータは、Julierの古典的な推奨値 \(\kappa = 3-n\) （\(n=2\) のとき \(\kappa=1\) ）と \(\alpha=1\) を使います（元記事にあった \(\alpha=0.5,\kappa=0\) という組み合わせがなぜ問題なのかは、後述の「落とし穴」で検証します）。

alpha = 1.0
kappa = 1.0  # Julierの古典的な推奨値 kappa = 3 - n（n=2 のとき kappa=1）

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

# シグマ点を格納する配列 (2n+1個のシグマ点)
sigma_points = np.zeros((n, 2 * n + 1))
# 式1: 最初のシグマ点は平均自身
sigma_points[:, 0] = x_mean
# (n+lambda) * Px の平方根（固有値分解ベースの sym_sqrt を使用）
sqrt_term = sym_sqrt((n + lam) * x_P)
# 式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]

for i in range(2 * n + 1):
    print(f"  sigma_{i}: {sigma_points[:, i]}")

# 重み関数を格納する配列
w = np.zeros(2 * n + 1)
# 式6, 7: 重み関数の計算
w[0] = lam / (n + lam)
w[1:] = 1 / (2 * (n + lam))
print(f"\nweights: {w}, sum: {w.sum()}")

# 変換されたシグマ点を格納する配列
sigma_y = np.zeros((m, 2 * n + 1))
# 式5: 非線形変換
for i in range(2 * n + 1):
    sigma_y[:, i] = f(sigma_points[:, i])

# 式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}")

実行結果:

lambda: 1.0, n+lambda: 3.0
  sigma_0: [0. 1.]
  sigma_1: [0.774597 2.549193]
  sigma_2: [1.549193 4.098387]
  sigma_3: [-0.774597 -0.549193]
  sigma_4: [-1.549193 -2.098387]

weights: [0.33333333 0.16666667 0.16666667 0.16666667 0.16666667], sum: 1.0

出力yの推定平均: [2.]
出力yの推定共分散行列:
[[5.16]]

証明(a)〜(c)が実際に成り立っていることを、シグマ点から入力の平均・共分散を再構成して確認します。

recon_mean = np.sum(w * sigma_points, axis=1)
recon_P = np.zeros((n, n))
for i in range(2 * n + 1):
    diff = sigma_points[:, i] - recon_mean
    recon_P += w[i] * np.outer(diff, diff)

print(f"再構成された平均: {recon_mean} (一致: {np.allclose(recon_mean, x_mean)})")
print(f"再構成された共分散:\n{recon_P}\n(一致: {np.allclose(recon_P, x_P)})")

実行結果:

再構成された平均: [0. 1.] (一致: True)
再構成された共分散:
[[1. 2.]
 [2. 4.]]
(一致: True)

証明の通り、シグマ点から再構成した平均・共分散が入力の \(\bar x\) , \(P_x\) と機械精度で一致しています。

4. 数値検証：解析解・モンテカルロ法との比較

\(X=(X_1,X_2)\) がガウス分布 \(\mathcal N(\bar x, P_x)\) に従う場合、\(Y=X_1X_2\) の平均・分散はウィックの定理（Isserlisの定理）を使って解析的に計算できます。\(X_1=\mu_1+\varepsilon_1\) , \(X_2=\mu_2+\varepsilon_2\) （\(\varepsilon\) は平均0のガウス変数、共分散を \(\Sigma\) とする）とおくと、

\[ E[Y] = \mu_1\mu_2 + \Sigma_{12} \] \[ \mathrm{Var}[Y] = \mu_1^2\Sigma_{22} + \mu_2^2\Sigma_{11} + \Sigma_{11}\Sigma_{22} + \Sigma_{12}^2 + 2\mu_1\mu_2\Sigma_{12} \]

（分散の式は、4次モーメントに関するIsserlisの定理 \(E[\varepsilon_1^2\varepsilon_2^2]=\Sigma_{11}\Sigma_{22}+2\Sigma_{12}^2\) と、ガウス変数の奇数次中心モーメントが恒等的に0になる性質から導出されます）。本記事の数値例に代入して検証します。

mu1, mu2 = x_mean
S11, S22, S12 = x_P[0, 0], x_P[1, 1], x_P[0, 1]
true_mean = mu1 * mu2 + S12
true_var = mu1**2 * S22 + mu2**2 * S11 + S11 * S22 + S12**2 + 2 * mu1 * mu2 * S12
print(f"解析解: true_mean = {true_mean}, true_var = {true_var}")

rng = np.random.default_rng(0)
N = 2_000_000
# x_Pは特異行列なので method="eigh"（固有値分解ベースのサンプリング）を使う
samples = rng.multivariate_normal(x_mean, x_P, size=N, method="eigh")
y_samples = samples[:, 0] * samples[:, 1]
print(f"モンテカルロ (N={N}): 平均={y_samples.mean():.4f}, 分散={y_samples.var():.4f}")

実行結果:

解析解: true_mean = 2.0, true_var = 9.0
モンテカルロ (N=2000000): 平均=2.0010, 分散=9.0168
手法平均 \(E[Y]\)分散 \(\mathrm{Var}[Y]\)
UT（\(\alpha=1,\kappa=1\) ）2.05.16
解析解（ガウス仮定）2.09.0
モンテカルロ法（\(N=2\times10^6\) ）2.00109.0168
EKF相当（線形化のみ、\(f(\bar x)\) ）0.0―

証明の通り、UTの平均推定値2.0は、EKF相当の0.0とは対照的に、解析解・モンテカルロ法と完全に一致しています（\(f\) が2次関数であるため、証明で示した3次以降の残余項 \(R_3\) が恒等的に0になるからです）。一方、UTの分散推定値5.16は、真の値9.0を43%ほど過小評価しています。これは、\(\mathrm{Var}[Y]\) が \(X\) の4次モーメント（\(\Sigma_{11}\Sigma_{22}+\Sigma_{12}^2\) の項）に依存する量であり、証明(a)〜(c)が保証するのはあくまで入力の2次モーメント（平均・共分散）までの厳密な再現だからです。\(2n+1\) 個のシグマ点だけでは一般に4次モーメントまで再現できないため、分散の推定精度はスケーリングパラメータ \(\lambda\) （ひいては \(\alpha,\kappa\) ）の選び方に依存する近似にとどまります。

分布の形に依存しないのは平均だけであることの確認

証明で述べたとおり、2次関数 \(f\) に対する平均の厳密性は分布の形（ガウス分布かどうか）によらず成り立つはずです。これを、\(X_1,X_2\) が同じ平均・共分散を持ちながらガウス分布ではない、退化した2点分布で確認します（\(P_x\) が特異行列なので、非零固有値の固有ベクトル方向に2点を置くだけで同じ平均・共分散を再現できます）。

eigvals, eigvecs = np.linalg.eigh(x_P)
idx = np.argmax(eigvals)
direction = eigvecs[:, idx]
scale = np.sqrt(eigvals[idx])

pts = np.array([x_mean + scale * direction, x_mean - scale * direction])
probs = np.array([0.5, 0.5])

disc_mean = (probs[:, None] * pts).sum(axis=0)
diffs = pts - disc_mean
disc_cov = sum(p * np.outer(d, d) for p, d in zip(probs, diffs))
y_disc = pts[:, 0] * pts[:, 1]
disc_y_mean = (probs * y_disc).sum()
disc_y_var = (probs * (y_disc - disc_y_mean) ** 2).sum()

print(f"2点分布の平均: {disc_mean} (一致: {np.allclose(disc_mean, x_mean)})")
print(f"2点分布の共分散:\n{disc_cov}\n(一致: {np.allclose(disc_cov, x_P)})")
print(f"2点分布での E[Y]={disc_y_mean}, Var[Y]={disc_y_var}")

実行結果:

2点分布の平均: [0. 1.] (一致: True)
2点分布の共分散:
[[1. 2.]
 [2. 4.]]
(一致: True)
2点分布での E[Y]=2.0, Var[Y]=1.0

平均・共分散が全く同じ（\(\bar x=[0,1]\) , \(P_x=\begin{pmatrix}1&2\\2&4\end{pmatrix}\) ）にもかかわらず、この2点分布での \(\mathrm{Var}[Y]\) は1.0であり、ガウス分布を仮定した場合の9.0とは全く異なります。一方 \(E[Y]\) はどちらも2.0で完全に一致します。これは証明の帰結そのものです——2次関数 \(f\) の期待値は入力の平均・共分散だけで一意に定まりますが、分散（4次モーメントを含む量）は分布の形に依存するため、UTがどれだけ精密にシグマ点を配置しても、入力分布の4次以上のモーメント情報なしに厳密な分散を得ることはできません。

5. 図：シグマ点の配置と非線形変換後の分布

先ほどの例は \(P_x\) が特異行列（相関係数1）だったため、2次元平面上では1本の直線に縮退しており可視化には向きません。ここでは非特異な共分散行列を持つ別の例で、シグマ点の配置と非線形変換後の分布を可視化します。

x_P2 = np.array([[1.0, 1.5], [1.5, 4.0]])
print(f"det(P_x2) = {np.linalg.det(x_P2):.4f}, eigvals = {np.linalg.eigvalsh(x_P2)}")

sqrt_term2 = sym_sqrt((n + lam) * x_P2)
sigma_points2 = np.zeros((n, 2 * n + 1))
sigma_points2[:, 0] = x_mean
for i in range(n):
    sigma_points2[:, i + 1] = x_mean + sqrt_term2[:, i]
    sigma_points2[:, i + n + 1] = x_mean - sqrt_term2[:, i]

sigma_y2 = np.zeros((m, 2 * n + 1))
for i in range(2 * n + 1):
    sigma_y2[:, i] = f(sigma_points2[:, i])

y_mean2 = np.sum(w * sigma_y2, axis=1)
y_P2 = np.zeros((m, m))
for i in range(2 * n + 1):
    diff = sigma_y2[:, i] - y_mean2
    y_P2 += w[i] * np.outer(diff, diff)

S11b, S22b, S12b = x_P2[0, 0], x_P2[1, 1], x_P2[0, 1]
true_mean2 = mu1 * mu2 + S12b
true_var2 = mu1**2 * S22b + mu2**2 * S11b + S11b * S22b + S12b**2 + 2 * mu1 * mu2 * S12b

samples2 = rng.multivariate_normal(x_mean, x_P2, size=N)
y_samples2 = samples2[:, 0] * samples2[:, 1]

print(f"UT: y_mean2={y_mean2}, y_P2={y_P2}")
print(f"解析解: true_mean2={true_mean2}, true_var2={true_var2}")
print(f"モンテカルロ: 平均={y_samples2.mean():.4f}, 分散={y_samples2.var():.4f}")

実行結果:

det(P_x2) = 1.7500, eigvals = [0.378682 4.621318]
UT: y_mean2=[1.5], y_P2=[[2.644608]]
解析解: true_mean2=1.5, true_var2=7.25
モンテカルロ: 平均=1.5015, 分散=7.2674

この例でも、UTの平均推定値1.5は解析解・モンテカルロ法と完全に一致する一方、分散推定値2.64は真の値7.25を大きく下回っています。以下の図は、この非特異な例における5個のシグマ点の配置（左）と、非線形変換後の出力分布（右）を可視化したものです（可視化専用に \(N=20{,}000\) のモンテカルロサンプルを別途生成しています）。

シグマ点の配置と非線形変換後の分布：UTの平均は厳密に一致し分散は過小評価される

左図では、たった5個のシグマ点（青丸・青ひし形）が、20,000個中1,500個を描画したモンテカルロサンプル（緑の点群）が張る95%共分散楕円（灰色破線）を過不足なく代表していることがわかります。右図では、変換後の出力 \(Y=X_1X_2\) のヒストグラム（緑）に対して、UTの平均推定値（青破線）が真の平均（オレンジ実線）と完全に重なる一方、UTの \(\pm1\sigma\) 幅（青帯）が真の \(\pm1\sigma\) 幅（オレンジ帯）よりも明らかに狭く、分散を過小評価している様子が視覚的に確認できます。

落とし穴：負の重みと不正な共分散推定

式(6)を見ると、\(\lambda<0\) のとき中心シグマ点の重み \(w_0\) は負になり得ます。実は、元記事で使われていた \(\alpha=0.5,\kappa=0\) という組み合わせは、この落とし穴を踏んでいます。

alpha_bad, kappa_bad = 0.5, 0.0
lam_bad = alpha_bad**2 * (n + kappa_bad) - n
w_bad = np.zeros(2 * n + 1)
w_bad[0] = lam_bad / (n + lam_bad)
w_bad[1:] = 1 / (2 * (n + lam_bad))

sqrt_term_bad = sym_sqrt((n + lam_bad) * x_P)
sigma_points_bad = np.zeros((n, 2 * n + 1))
sigma_points_bad[:, 0] = x_mean
for i in range(n):
    sigma_points_bad[:, i + 1] = x_mean + sqrt_term_bad[:, i]
    sigma_points_bad[:, i + n + 1] = x_mean - sqrt_term_bad[:, i]

sigma_y_bad = np.zeros((m, 2 * n + 1))
for i in range(2 * n + 1):
    sigma_y_bad[:, i] = f(sigma_points_bad[:, i])

y_mean_bad = np.sum(w_bad * sigma_y_bad, axis=1)
y_P_bad = np.zeros((m, m))
for i in range(2 * n + 1):
    diff = sigma_y_bad[:, i] - y_mean_bad
    y_P_bad += w_bad[i] * np.outer(diff, diff)

print(f"lambda_bad: {lam_bad}, n+lambda_bad: {n + lam_bad}")
print(f"weights_bad: {w_bad}, sum: {w_bad.sum()}")
print(f"y_mean_bad: {y_mean_bad}, y_P_bad: {y_P_bad}")

実行結果:

lambda_bad: -1.5, n+lambda_bad: 0.5
weights_bad: [-3.  1.  1.  1.  1.], sum: 1.0
y_mean_bad: [2.], y_P_bad: [[-1.64]]

平均推定値は証明の通り2.0と正しく求まりますが、共分散推定値が -1.64 という負の値になってしまいました。分散が負になることは数学的にあり得ないため、これは明らかに不正な推定です。原因は、共分散が階数1の外積 \((y_{\sigma,i}-\bar y)(y_{\sigma,i}-\bar y)^T\) の重み付き和（式9）として計算される点にあります。中心シグマ点の重み \(w_0=-3\) が大きく負であるため、他の4点からの正の寄与を打ち消して余りある大きさになり、結果全体が負に転じてしまったのです。

この問題は次元 \(n\) が大きいほど深刻になり、 Cubature Kalman Filter（CKF）の理論とPython実装 で詳しく証明しているように、Julierの古典的な推奨則 \(\kappa=3-n\) を使うUKFでは \(n>9\) で \(w_0\) が必ず負になります。実務上の対策としては、

  • \(\alpha\) を1に近い値にする（\(\lambda\) が過度に負にならないようにする）、または
  • Unscented Kalman Filter（UKF）の理論とPython実装 で使われている van der Merwe のスケーリング（\(w_0^{(c)}=\lambda/(n+\lambda)+(1-\alpha^2+\beta)\) という平均用と共分散用で異なる重みを使う拡張）を採用する、

といった方法が挙げられます。いずれにせよ、UTの重みは常に検証すべきものであり、盲目的に式(6)(7)を適用すると不正な共分散が得られる場合があるという点は、実装上の重要な注意点です。

最新研究動向

UTの重みが分布の形状（特にガウス性）を暗黙に仮定した固定式であるという、まさに上記の落とし穴に関連する研究が近年も進んでいます。Majewski, Modzelewski, Żugaj, & Lichota (2026) の “Robust Unscented Kalman Filtering via Recurrent Meta-Adaptation of Sigma-Point Weights”（ arXiv:2603.04360 ）は、固定のスケーリングパラメータに基づく従来のシグマ点重み付けが非ガウス性や時変ダイナミクスに追従できないという課題に対し、シグマ点の重み付けそのものをハイパーパラメータ最適化問題として再定式化する Meta-Adaptive UKF（MA-UKF）を提案しています。観測イノベーションの履歴を再帰的なコンテキストエンコーダで圧縮し、その埋め込みからポリシーネットワークが各時刻の重みを動的に生成することで、重い裾を持つ観測ノイズや未知の動的レジームへの頑健性を高めています。これは本記事で示した「固定式の重みでは分散推定が分布依存で不正確になりうる」という限界に、データ駆動で対処しようとする試みだといえます。

また、UTは目標追跡・カルマンフィルタの枠を超えて、不確実性伝播の汎用ツールとしても引き続き活用されています。Chu, Shrestha, Gu, & Gross (2025) の “Robust Flower Cluster Matching Using The Unscented Transform”（ arXiv:2503.20631 ）は、農業用ロボットが果樹の花のクラスタをRGB-Dカメラで追跡する際、記述子の空間的な不確実性の許容範囲をUTで効率的に推定し、時間経過による見え方の変化や遮蔽に頑健な画像レジストレーションを実現しています。これは、本記事で証明した「非線形変換の前後で平均・共分散を少数の点だけで正確に伝播できる」というUTの基本性質が、2025年になってもロボティクスや農業分野の実応用で活用され続けていることを示す一例です。

本記事の位置づけ（関連記事との住み分け）

Unscented変換・シグマ点を扱う記事は複数ありますが、以下のように役割を分担しています。

  • 信号処理におけるフィルタリング手法の基礎 — スカラー変数を例に、UKFがテイラー展開の2次項まで捉えられる理由を直感的に説明するにとどめ、多変数での厳密な証明は本記事に委ねています。
  • 本記事 — 一般の次元 \(n\) において、シグマ点と重みの組が入力の平均・共分散（2次モーメント）を厳密に再現することを証明し、非線形変換後の平均が2次まで正確な理由をテイラー展開で導出。実行検証・図表・落とし穴・最新研究動向を扱います。
  • Cubature Kalman Filter（CKF）の理論とPython実装 — 球面キュバチャ則が3次の多項式まで厳密であることを証明しています。これは本記事の「2次モーメントまでの厳密性」とは異なる次数の主張であり、UKFとCKFの重みの安定性の違い（負の重み問題）についても定量的に検証しています。
  • Unscented Kalman Filter（UKF）の理論とPython実装 — 本記事とCKF記事の証明を前提として、カルマンフィルタへの実装（予測・更新ステップ）とEKFとの精度比較に焦点を当てています。

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

関連記事

関連ツール

この記事に関連する無料ツール

DevToolBox Moving Average Calculator SMA/EMA/WMA をブラウザで即座に計算 DevToolBox JSON Formatter JSONを整形・圧縮・バリデーション DevToolBox Regex Tester 正規表現をリアルタイムテスト

関連記事