💨

独立でない複数の確率変数の混合三角モーメントの計算(MATLABで確認)

2023/05/01に公開

はじめに

Moment-Based Kalman Filter(以降 MKF)というものが紹介されており，理解がてら試してみる．
詳細は下記の論文を参照．
MKFでは独立でない複数変数の混合三角モーメントの計算が必要となるが，
計算を進めるうえでいろいろな箇所で手こずったので，本記事では備忘がてら計算過程と結果を記す．
（やっていることは変数を直交基底に変数変換して，独立した変数にしているだけだと思うのですが...その計算が大変...）
変数の取り方や記号はなるべく引用文献に合わせる．

流れてきたtweet↓

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]$ の計算は下記参照．

例

下記の分布の確率変数に対して， $\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を求める
% このときのlambdaは固有値が対角成分となる対角行列なので，要素の取出はインデックスに注意
[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となる．

% 計算する　答：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

ログインするとコメントできます