Zenn
📝

混合三角多項式の確率モーメント(mixed-trigonometric-polynomial moment)計算とMATLABでの確認

2023/04/01に公開
math
statistics
matlab
robotics
確率モーメント
tech

はじめに

確率モーメントベースで状態の不確実性を伝播させるという手法が下記論文(Moment-Based Exact Uncertainty Propagation Through NOnlinear Stochastic Autonomous Systems)で紹介されており，手計算とモンテカルロ法での確認を行ったので紹介する．
https://arxiv.org/abs/2101.12490

これにより，ある時刻でのロボットの状態（位置など）から数ステップ後の状態を不確実性も含めて扱うことができる．
論文においては，混合三角多項式のモーメント計算が示されている．

本記事では，ある確率分布\thetaに従う確率モーメントを計算した結果を示す．
分布としては正規分布，一様分布，指数分布を扱う．

※原著の英語での用語を日本語にしているので，翻訳違いは容赦ください．

特性関数と確率モーメント

一般式

確率変数\thetaの特性関数を\Phi_\thetaとする．
混合三角多項式の\alpha次モーメント（\alpha=\alpha_1+\alpha_2+\alpha_3）は一般に下記式で計算できる．
(論文Lemma4)

\begin{align*} &\mathbb{E}[\theta^{\alpha_1}\cos^{\alpha_2}(\theta)\sin^{\alpha_3}(\theta)] \\ &=\frac{1}{i^{\alpha_1+\alpha_3}2^{\alpha_2+\alpha_3}}\sum_{(k_1,k_2)=(0,0)}^{(\alpha_2,\alpha_3)} {}_{\alpha_2}\mathrm{C}_{k_1}{}_{\alpha_3}\mathrm{C}_{k_2} (-1)^{\alpha_3-k_2} \left . \frac{\mathrm{d}^{\alpha_1}}{\mathrm{d}t^{\alpha_1}}\Phi_\theta(t) \right|_{t=2(k_1+k_2)-\alpha_2-\alpha_3} \end{align*}

混合三角多項式のモーメントの例

モーメントの例\Phi(t)(=\Phi_\theta(t))を用いて計算する．
{\Phi}'(t) および {\Phi}''(t) をそれぞれ一階微分，二階微分とすると以下のようになる．

\begin{align*} &\mathbb{E}[\cos(\theta)] &=& \frac{1}{2}\{\Phi(1)+\Phi(-1)\} \\ &\mathbb{E}[\sin(\theta)] &=& \frac{1}{2i}\{\Phi(1)-\Phi(-1)\} \\ &\mathbb{E}[\theta\cos(\theta)] &=& \frac{1}{2i}\{\Phi'(1)+\Phi'(-1)\} \\ &\mathbb{E}[\theta\sin(\theta)] &=& -\frac{1}{2}\{\Phi'(1)-\Phi'(-1)\} \\ &\mathbb{E}[\theta^2\cos(\theta)] &=& -\frac{1}{2}\{\Phi''(1)+\Phi''(-1)\} \\ &\mathbb{E}[\theta^2\sin(\theta)] &=& -\frac{1}{2i}\{\Phi''(1)-\Phi''(-1)\} \end{align*}

様々な分布でのモーメントの計算

正規分布，一様分布，指数分布の確率変数\thetaについてそれぞれ計算する．
また，モンテカルロ法により計算結果を確かめる．

モンテカルロ法による計算方法

モンテカルロ法では，確率分布に従う乱数\thetaをN個（10^8個）取り出し，計算結果の平均をとる．
例えば，\rm{average}(\cos(\theta)) \approx \mathbb{E}[\cos(\theta)]とする

各確率分布に従う乱数は下記のように取得できる（MATLAB）．
参考：乱数 - MATLAB random- MathWorks 日本

%　乱数のシード値の固定
rng('default');
rng(1);
% 正規分布の乱数(平均μ=0.5，標準偏差σ=0.4)
rng(1);
theta = random('Normal',0.5,0.4,[1,1e8]);
% 一様分布の乱数(下限-0.2,上限0.5)
rng(1);
theta = random('Uniform',-0.2,0.5,[1,1e8]);
% 指数分布の乱数(平均値1/3)
rng(1);
theta = random('Exponential',1/3.0,[1,1e8]);

次に，確率モーメントを計算する．

% E[cos(θ)]
mean(cos(theta))
% E[sin(θ)]
mean(sin(theta))
% E[θcos(θ)]
mean(theta.*cos(theta))
% E[θsin(θ)]
mean(theta.*sin(theta))
% E[θ^2 cos(θ)]
mean(theta.*theta.*cos(theta))
% E[θ^2 sin(θ)]
mean(theta.*theta.*sin(theta))

以上により各分布でのモーメントが近似できる．

特性関数からの計算と確認

確率モーメントについて手計算する．
またモンテカルロ法での近似解との比較を行う

正規分布

平均値\mu，分散\sigma^2とする．
特性関数は

\Phi_\theta(t)=\exp\left(i\mu t-\frac{\sigma^2 t^2}{2}\right)

である．

\begin{align*} &\mathbb{E}[\cos(\theta)] &=& \cos(\mu)\exp\left(-{\sigma^2}/{2}\right)\\ &\mathbb{E}[\sin(\theta)] &=& \sin(\mu)\exp\left(-{\sigma^2}/{2}\right) \\ &\mathbb{E}[\theta\cos(\theta)] &=& (\mu \cos (\mu)-\sigma^2 \sin (\mu))\exp\left(-{\sigma^2}/{2}\right) \\ &\mathbb{E}[\theta\sin(\theta)] &=& (\mu \sin (\mu)+\sigma^2 \cos (\mu))\exp\left(-{\sigma^2}/{2}\right) \\ &\mathbb{E}[\theta^2\cos(\theta)] &=& \{(\sigma^2+\mu^2-\sigma^4)\cos(\mu)-2\mu\sigma^2\sin(\mu)\} \exp\left(-{\sigma^2}/{2}\right) \\ &\mathbb{E}[\theta^2\sin(\theta)] &=& \{(\sigma^2+\mu^2-\sigma^4)\sin(\mu)+2\mu\sigma^2\cos(\mu)\}\exp\left(-{\sigma^2}/{2}\right) \end{align*}

モンテカルロ法との比較は下記となり，計算が合っていることが確かめられる．

\mathbb{E}[\cos(\theta)] \mathbb{E}[\sin(\theta)] \mathbb{E}[\theta\cos(\theta)] \mathbb{E}[\theta\sin(\theta)] \mathbb{E}[\theta^2\cos(\theta)] \mathbb{E}[\theta^2\sin(\theta)]
解析解 0.8101 0.4426 0.3342 0.3509 0.2406 0.2997
モンテカルロ 0.8101 0.4426 0.3343 0.3509 0.2406 0.2998
差( \times 10^{-5} ) 1.3051 -2.3970 -1.4950 -2.4924 -1.9404 -2.7647
コード
mu = 0.5;
sigma = 0.4;
Ecos = cos(mu)*exp(-sigma^2/2)
Esin = sin(mu)*exp(-sigma^2/2)
Ethcos = (mu*cos(mu)-sigma^2*sin(mu))*exp(-sigma^2/2)
Ethsin = (mu*sin(mu)+sigma^2*cos(mu))*exp(-sigma^2/2)
Eth2cos = ((sigma^2+mu^2-sigma^4)*cos(mu)-2*mu*sigma^2*sin(mu))*exp(-sigma^2/2)
Eth2sin = ((sigma^2+mu^2-sigma^4)*sin(mu)+2*mu*sigma^2*cos(mu))*exp(-sigma^2/2)

%% Monte Carlo   
rng('default');
rng(1);
theta = random('Normal',mu,sigma,[1,1e8]);

%% ====ここから以降は共通===
% E[cos(θ)]
mean(cos(theta))
% E[sin(θ)]
mean(sin(theta))
% E[θcos(θ)]
mean(theta.*cos(theta))
% E[θsin(θ)]
mean(theta.*sin(theta))
% E[θ^2 cos(θ)]
mean(theta.*theta.*cos(theta))
% E[θ^2 sin(θ)]
mean(theta.*theta.*sin(theta))

%%error
e1 = Ecos - mean(cos(theta))
e2 = Esin - mean(sin(theta))
e3 = Ethcos - mean(theta.*cos(theta))
e4 = Ethsin - mean(theta.*sin(theta))
e5 = Eth2cos - mean(theta.*theta.*cos(theta))
e6 = Eth2sin - mean(theta.*theta.*sin(theta))

%表示
s1 = sprintf('| | $\\\\mathbb{E}[\\\\cos(\\\\theta)]$ | $\\\\mathbb{E}[\\\\sin(\\\\theta)]$ | $\\\\mathbb{E}[\\\\theta\\\\cos(\\\\theta)]$ | $\\\\mathbb{E}[\\\\theta\\\\sin(\\\\theta)]$ | $\\\\mathbb{E}[\\\\theta^2\\\\cos(\\\\theta)]$ | $\\\\mathbb{E}[\\\\theta^2\\\\sin(\\\\theta)]$ |\n');
s2 = sprintf('| ---- | ---- | ---- | ---- | ---- | ---- | ---- |\n');
s3 = sprintf('|解析解|%.4f|%.4f|%.4f|%.4f|%.4f|%.4f|\n',Ecos,Esin,Ethcos,Ethsin,Eth2cos,Eth2sin);
s4 = sprintf('|モンテカルロ|%.4f|%.4f|%.4f|%.4f|%.4f|%.4f|\n',mean(cos(theta)),mean(sin(theta)),mean(theta.*cos(theta)),mean(theta.*sin(theta)),mean(theta.*theta.*cos(theta)),mean(theta.*theta.*sin(theta)));
s5 = sprintf('|差( $\\\\times 10^{-5}$ )|%.4f|%.4f|%.4f|%.4f|%.4f|%.4f|\n',e1*1e5,e2*1e5,e3*1e5,e4*1e5,e5*1e5,e6*1e5);
fprintf([s1,s2,s3,s4,s5])
一様分布

分布の下限をL，上限をUとする．
特性関数は

\Phi_\theta(t)=\frac{\exp(iUt)-\exp(iLt)}{i(U-L)t}

である．

\begin{align*} &\mathbb{E}[\cos(\theta)] &=& \frac{\sin(U)-\sin(L)}{U-L}\\ &\mathbb{E}[\sin(\theta)] &=& -\frac{\cos(U)-\cos(L)}{U-L}\\ &\mathbb{E}[\theta\cos(\theta)] &=& -\frac{-U\sin(U)+L\sin(L)-\cos(U)+\cos(L)}{U-L}\\ &\mathbb{E}[\theta\sin(\theta)] &=& -\frac{U\cos(U)-L\cos(L)-\sin(U)+\sin(L)}{U-L}\\ &\mathbb{E}[\theta^2\cos(\theta)] &=& -\frac{(-U^2+2)\sin(U)-2U\cos(U)+(L^2-2)\sin(L)+2L\cos(L)}{U-L}\\ &\mathbb{E}[\theta^2\sin(\theta)] &=& \frac{(-U^2+2)\cos(U)+2U\sin(U)+(L^2-2)\cos(L)-2L\sin(L)}{U-L} \end{align*}

モンテカルロ法との比較は下記．

\mathbb{E}[\cos(\theta)] \mathbb{E}[\sin(\theta)] \mathbb{E}[\theta\cos(\theta)] \mathbb{E}[\theta\sin(\theta)] \mathbb{E}[\theta^2\cos(\theta)] \mathbb{E}[\theta^2\sin(\theta)]
解析解 0.9687 0.1464 0.1393 0.0618 0.0589 0.0211
モンテカルロ 0.9687 0.1464 0.1393 0.0619 0.0589 0.0211
差( \times 10^{-5} ) 0.5153 -2.2585 -2.0843 -1.0109 -0.9325 -0.5151
コード
L = -0.2;
U = 0.5;

Ecos = (sin(U)-sin(L))/(U-L)
Esin = -(cos(U)-cos(L))/(U-L)
Ethcos = -(-U*sin(U)+L*sin(L)-cos(U)+cos(L))/(U-L)
Ethsin = -(U*cos(U)-L*cos(L)-sin(U)+sin(L))/(U-L)
Eth2cos =-((-U^2+2)*sin(U)-2*U*cos(U)+(L^2-2)*sin(L)+2*L*cos(L))/(U-L)
Eth2sin =((-U^2+2)*cos(U)+2*U*sin(U)+(L^2-2)*cos(L)-2*L*sin(L))/(U-L)

rng(1);
theta = random('Uniform',L,U,[1,1e8]);

%~~~以降省略~~~
指数分布

平均値を1/\lambdaとする．
特性関数は

\Phi_\theta(t)=\frac{\lambda}{\lambda-it}

である．

\begin{align*} &\mathbb{E}[\cos(\theta)] &=& \frac{\lambda^2}{\lambda^2+1}\\ &\mathbb{E}[\sin(\theta)] &=& \frac{\lambda}{\lambda^2+1}\\ &\mathbb{E}[\theta\cos(\theta)] &=& \frac{\lambda^3-\lambda}{(\lambda^2+1)^2}\\ &\mathbb{E}[\theta\sin(\theta)] &=& \frac{2\lambda^2}{(\lambda^2+1)^2}\\ &\mathbb{E}[\theta^2\cos(\theta)] &=& \frac{2\lambda^2(\lambda^2-3)}{(\lambda^2+1)^3}\\ &\mathbb{E}[\theta^2\sin(\theta)] &=& \frac{2\lambda(3\lambda^2-1)}{(\lambda^2+1)^3}\\ \end{align*}

モンテカルロ法との比較は下記．

\mathbb{E}[\cos(\theta)] \mathbb{E}[\sin(\theta)] \mathbb{E}[\theta\cos(\theta)] \mathbb{E}[\theta\sin(\theta)] \mathbb{E}[\theta^2\cos(\theta)] \mathbb{E}[\theta^2\sin(\theta)]
解析解 0.9000 0.3000 0.2400 0.1800 0.1080 0.1560
モンテカルロ 0.9000 0.3000 0.2400 0.1800 0.1080 0.1560
差( \times 10^{-5} ) -1.4157 2.2196 0.7986 2.0447 -0.7383 3.0397
コード
lambda = 3;

Ecos = lambda^2/(lambda^2+1)
Esin = lambda/(lambda^2+1)
Ethcos = (lambda^3-lambda)/(lambda^2+1)^2
Ethsin = (2*lambda^2)/(lambda^2+1)^2
Eth2cos =(2*lambda^2)*(lambda^2-3)/(lambda^2+1)^3
Eth2sin =(2*lambda)*(3*lambda^2-1)/(lambda^2+1)^3

rng(1);
theta = random('Exponential',1/lambda,[1,1e8]);

%~~~以降省略~~~

おわりに

モーメントの手計算を行ってみたが，かなり大変だという印象を受けたので，もし活用していただけるのであれば幸いである．
（特性関数の二階微分とか大変…）

Discussion