はじめに
Moment-Based Kalman Filter(以降 MKF)というものが紹介されており，理解がてら試してみる．
詳細は下記の論文を参照．
MKFでは独立でない複数変数の混合三角モーメントの計算が必要となるが，
計算を進めるうえでいろいろな箇所で手こずったので，本記事では備忘がてら計算過程と結果を記す．
（やっていることは変数を直交基底に変数変換して，独立した変数にしているだけだと思うのですが...その計算が大変...）
変数の取り方や記号はなるべく引用文献に合わせる．
流れてきたtweet↓
https://twitter.com/purewater0901/status/1617706889642639360?s=20
MKFについての論文↓
https://arxiv.org/abs/2301.09130
混合三角モーメント
本記事で扱う三角モーメントとは，例えば下記のようなものとなる．
ここで，xと\thetaは互いに独立でない確率変数である．
\begin{align*}
&\mathbb{E}[x\cos(\theta)]
\end{align*}
前提
通常の変数であれば上記の式は計算できないが，引用論文おいては，
確率変数が正規分布に従う場合には計算可能となる．
そこで，以降では各確率変数は正規分布に従うものとする．
\bm{x}=[x,\theta]^\mathrm{T}としたとき，
\begin{align*}
&\bm{x} \sim \mathcal{N} \left(\mu_{\bm{x}},\Sigma_{\bm{x}} \right)
\end{align*}
となる．
\mu_{\bm{x}}=[\mu_{{x}},\mu_{\theta}]^{\mathrm{T}}であり，xと\thetaの平均値を表し，\Sigma_{\bm{x}}は共分散行列である．
共分散行列については下記がわかりやすいと思う．
https://manabitimes.jp/math/1020
計算が大変なので2変数で表記しているが，もっと多くても基本的には同じ．
変数変換
はじめに，共分散行列の固有値\lambda_1,\lambda_2と直交行列Tを求める．Tの求め方は調べれば出てくるので割愛する．
このとき，
\begin{align*}
&T^{-1}\Sigma_{\bm{x}}T = \Lambda \\
&\Lambda = \begin{bmatrix}
\lambda_1&0\\
0&\lambda_2
\end{bmatrix}
\end{align*}
である．
次に，\bm{y}=T^{-1}\bm{x}とおくと，\mu_{\bm{y}}=T^{-1}\mu_{\bm{x}}=T^{\mathrm{T}}\mu_{\bm{x}}，\Sigma_{\bm{y}}=\Lambdaとなる．
ここの計算の詳細は引用文献の(5)式を参照．
このとき，下記のように表記される．
\begin{align*}
&\bm{y} \sim \mathcal{N} \left(\mu_{\bm{y}},\Sigma_{\bm{y}} \right)
\end{align*}
この変換を施すことにより，確率変数を独立した変数に変換できるため，混合三角モーメントを計算できるようになる．
モーメントの計算の一例
最初に挙げた\mathbb{E}[x\cos(\theta)]を対象に計算をしてみる．
引用文献のExample1と同様であるが計算過程は異なる手順をとる(なお，引用文献の計算は誤植と思われ，公開されているソースコードの中身は同じ計算結果となっている)．
先ほどの\bm{y}を下記のように表記しなおす．
\begin{align*}
&\bm{y}=(y_1,y_2)^{T} \sim \mathcal{N} \left(\mu_{\bm{y}},\Sigma_{\bm{y}} \right) = \mathcal{N}
\left(
\begin{pmatrix}
\bar{y_1}\\\bar{y_2}
\end{pmatrix},
\begin{pmatrix}
\lambda_1&0\\ 0&\lambda_2\\
\end{pmatrix}
\right)
\end{align*}
また,T_{ij}はTの(i,j)要素とする．
これを用いて計算すると下記のようになる．
\begin{align*}
\mathbb{E}&[x \cos(\theta)] \\
=& \mathbb{E}[(T_{11}y_1+T_{12}y_2)\cos(T_{21}y_1+T_{22}y_2)] \\
=&\mathbb{E}[(T_{11}y_1+T_{12}y_2)\{\cos(T_{21}y_1)\cos(T_{22}y_2)-\sin(T_{21}y_1)\sin(T_{22}y_2)\}] //加法定理 \\
=&\mathbb{E}[T_{11}y_1\cos(T_{21}y_1)\cos(T_{22}y_2)+T_{12}y_2\cos(T_{21}y_1)\cos(T_{22}y_2) \\
&-T_{11}y_1\sin(T_{21}y_1)\sin(T_{22}y_2)-(T_{12}y_2)\sin(T_{21}y_1)\sin(T_{22}y_2)] \\
=&\mathbb{E}
\left\lbrack
\frac{T_{11}}{T_{21}}\{T_{21}y_1\cos(T_{21}y_1)\cos(T_{22}y_2)\}+
\frac{T_{12}}{T_{22}}\{T_{22}y_2\cos(T_{21}y_1)\cos(T_{22}y_2)\} \right.\\
&\left. -\frac{T_{11}}{T_{21}}\{T_{21}y_1\sin(T_{21}y_1)\sin(T_{22}y_2)\}
-\frac{T_{12}}{T_{22}}\{T_{22}y_2\sin(T_{21}y_1)\sin(T_{22}y_2)\}
\right\rbrack \\
=&\mathbb{E}
\left\lbrack
\frac{T_{11}}{T_{21}}\{T_{21}y_1\cos(T_{21}y_1)\cos(T_{22}y_2)
-T_{21}y_1\sin(T_{21}y_1)\sin(T_{22}y_2)\}\right.\\
&\left. +\frac{T_{12}}{T_{22}}\{T_{22}y_2\cos(T_{21}y_1)\cos(T_{22}y_2)
-T_{22}y_2\sin(T_{21}y_1)\sin(T_{22}y_2)\}
\right\rbrack \\
\end{align*}
ここで文献と同様に，l_1=T_{21}y_1，l_2=T_{22}y_2とおくと
l_1 \sim \mathcal{N}(T_{21}\bar{y_1},T_{21}^2\lambda_1)およびl_2 \sim \mathcal{N}(T_{22}\bar{y_2},T_{22}^2\lambda_2)となっており，これらは独立した変数となる．
\begin{align*}
\mathbb{E}&[x \cos(\theta)] \\
=&\mathbb{E}
\left\lbrack
\frac{T_{11}}{T_{21}}\{l_1\cos(l_1)\cos(l_2)
-l_1\sin(l_1)\sin(l_2)\}\right.\\
&\left. +\frac{T_{12}}{T_{22}}\{l_2\cos(l_1)\cos(l_2)
-l_2\sin(l_1)\sin(l_2)\}
\right\rbrack \\
=&\frac{T_{11}}{T_{21}}
\{\mathbb{E}[l_1\cos(l_1)] \mathbb{E}[\cos(l_2)]
-\mathbb{E}[l_1\sin(l_1)]\mathbb{E}[\sin(l_2)]\} \\
&+\frac{T_{12}}{T_{22}}
\{\mathbb{E}[\cos(l_1)] \mathbb{E}[l_2\cos(l_2)]
-\mathbb{E}[\sin(l_1)]\mathbb{E}[l_2\sin(l_2)]\} \\
\end{align*}
上記の各\mathbb{E}[\cdot]は計算できるので，もともとの\mathbb{E}[x \cos(\theta)]を計算可能となる．
なお，各\mathbb{E}[\cdot]の計算は下記参照．
https://zenn.dev/spargel/articles/4a64d0157017cc
例
下記の分布の確率変数に対して，\mathbb{E}[x \cos(\theta)]を計算する（文献に合わせる）．
なお，対角化の計算などはMATLABを使用する（手計算面倒なので...）．
\begin{align*}
& \begin{pmatrix}
x\\ \theta
\end{pmatrix}
\sim \mathcal{N}
\left(
\begin{pmatrix}
{10.0}\\{\pi/3}
\end{pmatrix},
\begin{pmatrix}
5.0&1.5\\ 1.5&\pi/6\\
\end{pmatrix}
\right)
\end{align*}
まず，直交行列と固有値を求める．
mu_x = [10.0;pi/3]
Sigma_x = [5.0,1.5;1.5,pi/6]
[T,lambda] = eig(Sigma_x)
次に，うえで表記したように\bm{y}=T^{-1}\bm{x}とおいて，\mu_{\bm{y}}=T^{\mathrm{T}}\mu_{\bm{x}}とする．
また，l_1，l_2の平均値および分散も計算する．
mu_y = T'*mu_x;
mu_l1 = T(2,1)*mu_y(1);
mu_l2 = T(2,2)*mu_y(2);
sigma_l1 = T21^2*lambda(1,1)
sigma_l2 = T22^2*lambda(2,2)
そして，\mathbb{E}[l_1\cos(l_1)]などを計算する．
E_cosl1 = cos(mu_l1)*exp(-sigma_l1/2);
E_cosl2 = cos(mu_l2)*exp(-sigma_l2/2);
E_sinl1 = sin(mu_l1)*exp(-sigma_l1/2);
E_sinl2 = sin(mu_l2)*exp(-sigma_l2/2);
E_l1cosl1 = (mu_l1*cos(mu_l1)-sigma_l1*sin(mu_l1))*exp(-sigma_l1/2);
E_l2cosl2 = (mu_l2*cos(mu_l2)-sigma_l2*sin(mu_l2))*exp(-sigma_l2/2);
E_l1sinl1 = (mu_l1*sin(mu_l1)+sigma_l1*cos(mu_l1))*exp(-sigma_l1/2);
E_l2sinl2 = (mu_l2*sin(mu_l2)+sigma_l2*cos(mu_l2))*exp(-sigma_l2/2);
これを用いて\mathbb{E}[x \cos(\theta)]を計算すると，結果が2.8485となる．
T(1,1)/T(2,1)*(E_l1cosl1*E_cosl2-E_l1sinl1*E_sinl2)...
+T(1,2)/T(2,2)*(E_cosl1*E_l2cosl2-E_sinl1*E_l2sinl2)
上記の結果は引用文献と一致するので，計算は正しくできたといえるが，モンテカルロ法でも確認すると，2.8481となる．
当然シード値を変えると答えも変わるが，2.848までは同じなので，上記の結果は正しく計算できているといってよいかと思う．
rng('default');
rng(1);
x = mvnrnd(mu_x,Sigma_x,1e8);
mean(x(:,1).*cos(x(:,2)))
おわりに
本記事では混合三角モーメントの計算を書き下した．
記事がだいぶ長くなったので，ほかの形のモーメントの計算については次の記事に譲る．
Discussion