Zenn
📒

独立でない変数の混合三角モーメントの計算の色々な例

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

はじめに

これまでの記事（下記）で混合三角モーメントの計算を行ってきた．
そこで本記事では複数の確率変数の場合のモーメントについて計算する．
なお，混合三角モーメントはMixed-trigonometric-Polynominal Momentsの直訳であり本記事ではこの記載でする．

https://zenn.dev/spargel/articles/4a64d0157017cc

https://zenn.dev/spargel/articles/a2a5ce1b6cb421

また参考に引用している文献は下記．
https://arxiv.org/abs/2301.09130

MATLABで混合三角モーメント計算用クラスの作成

混合三角モーメントをそのたびに計算するのは面倒なので，MATLABのクラスにまとめる．
なお，対象とする確率分布は正規分布（ガウス分布）のみである．
クラス名はGaussianMTPMであり，引数は平均\muと\sigma^2ある．

MATLAB classコード
% compute E[f(x)],x~N(mu sigma2)
% input: mu,sigma (scalar) 
% include Mixed-trigonometric-Polynominal Moments
% E[\x^a1 cos^a2(x) sin^a3(x)],
% refference:https://arxiv.org/pdf/2101.12490.pdf
classdef GaussianMTPM
    properties
        mu
        sigma2
    end
    methods
        function obj = GaussianMTPM(mu,sigma2)
            if nargin == 2
                obj.mu = mu;
                obj.sigma2 = sigma2;
            end
        end
        function e = X(obj)
            e = obj.mu;
        end
        function e = X2(obj)
            % E[x^2]
            e = obj.sigma2+obj.mu^2;
        end
        function e = CosX(obj)
            % E[cosx]
            e = cos(obj.mu)*exp(-0.5*obj.sigma2);
        end
        function e = SinX(obj)
            % E[sinx]
            e = sin(obj.mu)*exp(-0.5*obj.sigma2);
        end
        function e = XCosX(obj)
            % E[xcosx]
            e = (obj.mu*cos(obj.mu)-obj.sigma2*sin(obj.mu)) * exp(-0.5*obj.sigma2);
        end
        function e = XSinX(obj)
            % E[xsinx]
            e = (obj.mu*sin(obj.mu)+obj.sigma2*cos(obj.mu)) * exp(-0.5*obj.sigma2);
        end 
        function e = X2CosX(obj)
            % E[x^2 cosx]
            e = ((obj.sigma2+obj.mu^2-obj.sigma2^2)*cos(obj.mu)-2*obj.mu*obj.sigma2*sin(obj.mu))*exp(-0.5*obj.sigma2);
        end 
        function e = X2SinX(obj)
            % E[x^2 sinx]
            e =((obj.sigma2+obj.mu^2-obj.sigma2^2)*sin(obj.mu)+2*obj.mu*obj.sigma2*cos(obj.mu))*exp(-0.5*obj.sigma2);
        end 
    end
end

使い方

w \sim \mathcal{N}(a,b)という正規分布（平均a,分散b=\sigma^2）の確率変数の宣言を下記のようにする．

w = GaussianMTPM(a,b)

例えば，\mathbb{E}[w\cos(w)]の場合，[\cdot]の中を呼び出すようにする．

w.XCosX

同様に\mathbb{E}[w^2\sin(w)]は

w.X2SinX

のように呼び出す．

色々な混合三角モーメント（2変数の場合）

\mathbb{x}=(x,\theta)^\mathrm{T}\sim \mathcal{N}\left(\mu_x,\Sigma_x \right)，について，\Sigma_xを対角化する直交行列をT，固有値を\lambda_1,\lambda_2とする．
y=T^{-1}xとすると，x=T_{11}y_1+T_{12}y_2,\theta=T_{21}y_1+T_{22}y_2となる．
また，\mathbb{y}\sim \mathcal{N}(\mu_y=T^{-1}\mu_x，\Sigma_y = \Lambda)である．
このとき，各期待値は下記のような形になる．

\begin{align*} \mathbb{E}&[x\theta] \\ &= T_{11}T_{21}\mathbb{E}[y_1^2]+(T_{11}T_{22}+T_{12}T_{21})\mathbb{E}[y_1]\mathbb{E}[y_2]+T_{12}T_{22}\mathbb{E}[y_2^2]\\ \mathbb{E}&[x\cos(\theta)] \\ &=\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)]\} \\ \mathbb{E}&[x\sin(\theta)] \\ &=\frac{T_{11}}{T_{21}} \{\mathbb{E}[l_1\sin(l_1)] \mathbb{E}[\cos(l_2)] +\mathbb{E}[l_1\cos(l_1)]\mathbb{E}[\sin(l_2)]\} \\ &+\frac{T_{12}}{T_{22}} \{\mathbb{E}[\sin(l_1)] \mathbb{E}[l_2\cos(l_2)] +\mathbb{E}[\cos(l_1)]\mathbb{E}[l_2\sin(l_2)]\} \\ \mathbb{E}&[x\cos(\theta)\sin(\theta)] \\ &=\mathbb{E}\left[ x\frac{\sin(2\theta)}{2} \right] =\frac{1}{2} \mathbb{E}\left[ x{\sin(2\theta)} \right]\\ &=\frac{T_{11}}{4T_{21}} \{\mathbb{E}[m_1\sin(m_1)] \mathbb{E}[\cos(m_2)] +\mathbb{E}[m_1\cos(m_1)]\mathbb{E}[\sin(m_2)]\} \\ &+\frac{T_{12}}{4T_{22}} \{\mathbb{E}[\sin(m_1)] \mathbb{E}[m_2\cos(m_2)] +\mathbb{E}[\cos(m_1)]\mathbb{E}[m_2\sin(m_2)]\} \\ \end{align*}

以前の記事と同様に，l_1=T_{21}y_1，l_2=T_{22}y_2とおいている．
また，m_1=2T_{21}y_1,m_2=2T_{22}y_2とおいて，
m_1\sim \mathcal{N}(2T_{21}\mu_{y_1},(2T_{21})^2\lambda_1)，m_2\sim \mathcal{N}(2T_{22}\mu_{y_2},(2T_{22})^2\lambda_2)である．
y,l,mについての期待値は前節の方法で計算できるので，この期待値を求めることが可能である．

色々な混合三角モーメント（3変数の場合）

３変数の場合も同様に計算できる．
x=(x,y,\theta)^\mathrm{T}であり，Tおよび\lambdaの定義は２変数の場合と同じである．
なおTは3\times3行列であり，固有値\lambda_iは３つになる．
また，z_{ij}=T_{ij}y_jとする．例えば，z_{23}=T_{23}y_3である．
このとき，

\begin{align*} z_{ij} \sim \mathcal{N}(T_{ij}\bar{y_j},T_{ij}^2\lambda_j) \end{align*}

である．
さらに，記述簡略化のため，C_{ij}=\cos(z_{ij}),S_{ij}=\sin(z_{ij})と記述する．

モーメントの計算結果は下記の通り．

\begin{align*} \mathbb{E}&[x\cos(\theta)]\\ &=\frac{T_{11}}{T_{31}} \{\mathbb{E}[z_{31}C_{31}]\mathbb{E}[C_{32}]\mathbb{E}[C_{33}]-\mathbb{E}[z_{31}C_{31}]\mathbb{E}[S_{32}]\mathbb{E}[S_{33}]-\mathbb{E}[z_{31}S_{31}]\mathbb{E}[C_{32}]\mathbb{E}[S_{33}]-\mathbb{E}[z_{31}S_{31}]\mathbb{E}[S_{32}]\mathbb{E}[C_{33}]\} \\ &~+\frac{T_{12}}{T_{32}} \{\mathbb{E}[C_{31}]\mathbb{E}[z_{32}C_{32}]\mathbb{E}[C_{33}]-\mathbb{E}[C_{31}]\mathbb{E}[z_{32}S_{32}]\mathbb{E}[S_{33}]-\mathbb{E}[S_{31}]\mathbb{E}[z_{32}C_{32}]\mathbb{E}[S_{33}]-\mathbb{E}[S_{31}]\mathbb{E}[z_{32}S_{32}]\mathbb{E}[C_{33}]\} \\ &~+\frac{T_{13}}{T_{33}} \{\mathbb{E}[C_{31}]\mathbb{E}[C_{32}]\mathbb{E}[z_{33}C_{33}]-\mathbb{E}[C_{31}]\mathbb{E}[S_{32}]\mathbb{E}[z_{33}S_{33}]-\mathbb{E}[S_{31}]\mathbb{E}[C_{32}]\mathbb{E}[z_{33}S_{33}]-\mathbb{E}[S_{31}]\mathbb{E}[S_{32}]\mathbb{E}[z_{33}C_{33}]\} \\ \mathbb{E}&[x\sin(\theta)]\\ &=\frac{T_{11}}{T_{31}} \{-\mathbb{E}[z_{31}S_{31}]\mathbb{E}[S_{32}]\mathbb{E}[S_{33}]+\mathbb{E}[z_{31}S_{31}]\mathbb{E}[C_{32}]\mathbb{E}[C_{33}]+\mathbb{E}[z_{31}C_{31}]\mathbb{E}[S_{32}]\mathbb{E}[C_{33}]+\mathbb{E}[z_{31}C_{31}]\mathbb{E}[C_{32}]\mathbb{E}[S_{33}]\} \\ &~+\frac{T_{12}}{T_{32}} \{-\mathbb{E}[S_{31}]\mathbb{E}[z_{32}S_{32}]\mathbb{E}[S_{33}]+\mathbb{E}[S_{31}]\mathbb{E}[z_{32}C_{32}]\mathbb{E}[C_{33}]+\mathbb{E}[C_{31}]\mathbb{E}[z_{32}S_{32}]\mathbb{E}[C_{33}]+\mathbb{E}[C_{31}]\mathbb{E}[z_{32}C_{32}]\mathbb{E}[S_{33}]\} \\ &~+\frac{T_{13}}{T_{33}} \{-\mathbb{E}[S_{31}]\mathbb{E}[S_{32}]\mathbb{E}[z_{33}S_{33}]+\mathbb{E}[S_{31}]\mathbb{E}[C_{32}]\mathbb{E}[z_{33}C_{33}]+\mathbb{E}[C_{31}]\mathbb{E}[S_{32}]\mathbb{E}[z_{33}C_{33}]+\mathbb{E}[C_{31}]\mathbb{E}[C_{32}]\mathbb{E}[z_{33}S_{33}]\} \\ \mathbb{E}&[x^2\cos(\theta)]\\ &=\left(\frac{T_{11}}{T_{31}}\right)^2 \{\mathbb{E}[z_{31}^2C_{31}]\mathbb{E}[C_{32}]\mathbb{E}[C_{33}]-\mathbb{E}[z_{31}^2C_{31}]\mathbb{E}[S_{32}]\mathbb{E}[S_{33}]-\mathbb{E}[z_{31}^2S_{31}]\mathbb{E}[C_{32}]\mathbb{E}[S_{33}]-\mathbb{E}[z_{31}^2S_{31}]\mathbb{E}[S_{32}]\mathbb{E}[C_{33}]\} \\ &~+\left(\frac{T_{12}}{T_{32}}\right)^2 \{\mathbb{E}[C_{31}]\mathbb{E}[z_{32}^2C_{32}]\mathbb{E}[C_{33}]-\mathbb{E}[C_{31}]\mathbb{E}[z_{32}^2S_{32}]\mathbb{E}[S_{33}]-\mathbb{E}[S_{31}]\mathbb{E}[z_{32}^2C_{32}]\mathbb{E}[S_{33}]-\mathbb{E}[S_{31}]\mathbb{E}[z_{32}^2S_{32}]\mathbb{E}[C_{33}]\} \\ &~+\left(\frac{T_{13}}{T_{33}}\right)^2 \{\mathbb{E}[C_{31}]\mathbb{E}[C_{32}]\mathbb{E}[z_{33}^2C_{33}]-\mathbb{E}[C_{31}]\mathbb{E}[S_{32}]\mathbb{E}[z_{33}^2S_{33}]-\mathbb{E}[S_{31}]\mathbb{E}[C_{32}]\mathbb{E}[z_{33}^2S_{33}]-\mathbb{E}[S_{31}]\mathbb{E}[S_{32}]\mathbb{E}[z_{33}^2C_{33}]\} \\ &~+2\frac{T_{11}}{T_{31}}\frac{T_{12}}{T_{32}} \{\mathbb{E}[z_{31}C_{31}]\mathbb{E}[z_{32}C_{32}]\mathbb{E}[C_{33}]-\mathbb{E}[z_{31}C_{31}]\mathbb{E}[z_{32}S_{32}]\mathbb{E}[S_{33}]-\mathbb{E}[z_{31}S_{31}]\mathbb{E}[z_{32}C_{32}]\mathbb{E}[S_{33}]-\mathbb{E}[z_{31}S_{31}]\mathbb{E}[z_{32}S_{32}]\mathbb{E}[C_{33}]\} \\ &~+2\frac{T_{12}}{T_{32}}\frac{T_{13}}{T_{33}} \{\mathbb{E}[C_{31}]\mathbb{E}[z_{32}C_{32}]\mathbb{E}[z_{33}C_{33}]-\mathbb{E}[C_{31}]\mathbb{E}[z_{32}S_{32}]\mathbb{E}[z_{33}S_{33}]-\mathbb{E}[S_{31}]\mathbb{E}[z_{32}C_{32}]\mathbb{E}[z_{33}S_{33}]-\mathbb{E}[S_{31}]\mathbb{E}[z_{32}S_{32}]\mathbb{E}[z_{33}C_{33}]\} \\ &~+2\frac{T_{13}}{T_{33}}\frac{T_{11}}{T_{31}} \{\mathbb{E}[z_{31}C_{31}]\mathbb{E}[C_{32}]\mathbb{E}[z_{33}C_{33}]-\mathbb{E}[z_{31}C_{31}]\mathbb{E}[S_{32}]\mathbb{E}[z_{33}S_{33}]-\mathbb{E}[z_{31}S_{31}]\mathbb{E}[C_{32}]\mathbb{E}[z_{33}S_{33}]-\mathbb{E}[z_{31}S_{31}]\mathbb{E}[S_{32}]\mathbb{E}[z_{33}C_{33}]\} \\ \mathbb{E}&[x^2\sin(\theta)]\\ &=\left(\frac{T_{11}}{T_{31}}\right)^2 \{-\mathbb{E}[z_{31}^2S_{31}]\mathbb{E}[S_{32}]\mathbb{E}[S_{33}]+\mathbb{E}[z_{31}^2S_{31}]\mathbb{E}[C_{32}]\mathbb{E}[C_{33}]+\mathbb{E}[z_{31}^2C_{31}]\mathbb{E}[S_{32}]\mathbb{E}[C_{33}]+\mathbb{E}[z_{31}^2C_{31}]\mathbb{E}[C_{32}]\mathbb{E}[S_{33}]\} \\ &~+\left(\frac{T_{12}}{T_{32}}\right)^2 \{-\mathbb{E}[S_{31}]\mathbb{E}[z_{32}^2S_{32}]\mathbb{E}[S_{33}]+\mathbb{E}[S_{31}]\mathbb{E}[z_{32}^2C_{32}]\mathbb{E}[C_{33}]+\mathbb{E}[C_{31}]\mathbb{E}[z_{32}^2S_{32}]\mathbb{E}[C_{33}]+\mathbb{E}[C_{31}]\mathbb{E}[z_{32}^2C_{32}]\mathbb{E}[S_{33}]\} \\ &~+\left(\frac{T_{13}}{T_{33}}\right)^2 \{-\mathbb{E}[S_{31}]\mathbb{E}[S_{32}]\mathbb{E}[z_{33}^2S_{33}]+\mathbb{E}[S_{31}]\mathbb{E}[C_{32}]\mathbb{E}[z_{33}^2C_{33}]+\mathbb{E}[C_{31}]\mathbb{E}[S_{32}]\mathbb{E}[z_{33}^2C_{33}]+\mathbb{E}[C_{31}]\mathbb{E}[C_{32}]\mathbb{E}[z_{33}^2S_{33}]\} \\ &~+2\frac{T_{11}}{T_{31}}\frac{T_{12}}{T_{32}} \{-\mathbb{E}[z_{31}S_{31}]\mathbb{E}[z_{32}S_{32}]\mathbb{E}[S_{33}]+\mathbb{E}[z_{31}S_{31}]\mathbb{E}[z_{32}C_{32}]\mathbb{E}[C_{33}]+\mathbb{E}[z_{31}C_{31}]\mathbb{E}[z_{32}S_{32}]\mathbb{E}[C_{33}]+\mathbb{E}[z_{31}C_{31}]\mathbb{E}[z_{32}C_{32}]\mathbb{E}[S_{33}]\} \\ &~+2\frac{T_{12}}{T_{32}}\frac{T_{13}}{T_{33}} \{-\mathbb{E}[S_{31}]\mathbb{E}[z_{32}S_{32}]\mathbb{E}[z_{33}S_{33}]+\mathbb{E}[S_{31}]\mathbb{E}[z_{32}C_{32}]\mathbb{E}[z_{33}C_{33}]+\mathbb{E}[C_{31}]\mathbb{E}[z_{32}S_{32}]\mathbb{E}[z_{33}C_{33}]+\mathbb{E}[C_{31}]\mathbb{E}[z_{32}C_{32}]\mathbb{E}[z_{33}S_{33}]\} \\ &~+2\frac{T_{13}}{T_{33}}\frac{T_{11}}{T_{31}} \{-\mathbb{E}[z_{31}S_{31}]\mathbb{E}[S_{32}]\mathbb{E}[z_{33}S_{33}]+\mathbb{E}[z_{31}S_{31}]\mathbb{E}[C_{32}]\mathbb{E}[z_{33}C_{33}]+\mathbb{E}[z_{31}C_{31}]\mathbb{E}[S_{32}]\mathbb{E}[z_{33}C_{33}]+\mathbb{E}[z_{31}C_{31}]\mathbb{E}[C_{32}]\mathbb{E}[z_{33}S_{33}]\} \\ \end{align*}

見た目は仰々しいが，三角関数の係数が同次数の場合は\{\}の前の係数が変わるだけなので，上記を計算してしまえば残りは比較的楽...

続き
\begin{align*} \mathbb{E}&[y\cos(\theta)]\\ &=\frac{T_{21}}{T_{31}} \{\mathbb{E}[z_{31}C_{31}]\mathbb{E}[C_{32}]\mathbb{E}[C_{33}]-\mathbb{E}[z_{31}C_{31}]\mathbb{E}[S_{32}]\mathbb{E}[S_{33}]-\mathbb{E}[z_{31}S_{31}]\mathbb{E}[C_{32}]\mathbb{E}[S_{33}]-\mathbb{E}[z_{31}S_{31}]\mathbb{E}[S_{32}]\mathbb{E}[C_{33}]\} \\ &~+\frac{T_{22}}{T_{32}} \{\mathbb{E}[C_{31}]\mathbb{E}[z_{32}C_{32}]\mathbb{E}[C_{33}]-\mathbb{E}[C_{31}]\mathbb{E}[z_{32}S_{32}]\mathbb{E}[S_{33}]-\mathbb{E}[S_{31}]\mathbb{E}[z_{32}C_{32}]\mathbb{E}[S_{33}]-\mathbb{E}[S_{31}]\mathbb{E}[z_{32}S_{32}]\mathbb{E}[C_{33}]\} \\ &~+\frac{T_{23}}{T_{33}} \{\mathbb{E}[C_{31}]\mathbb{E}[C_{32}]\mathbb{E}[z_{33}C_{33}]-\mathbb{E}[C_{31}]\mathbb{E}[S_{32}]\mathbb{E}[z_{33}S_{33}]-\mathbb{E}[S_{31}]\mathbb{E}[C_{32}]\mathbb{E}[z_{33}S_{33}]-\mathbb{E}[S_{31}]\mathbb{E}[S_{32}]\mathbb{E}[z_{33}C_{33}]\} \\ \mathbb{E}&[y\sin(\theta)]\\ &=\frac{T_{21}}{T_{31}} \{-\mathbb{E}[z_{31}S_{31}]\mathbb{E}[S_{32}]\mathbb{E}[S_{33}]+\mathbb{E}[z_{31}S_{31}]\mathbb{E}[C_{32}]\mathbb{E}[C_{33}]+\mathbb{E}[z_{31}C_{31}]\mathbb{E}[S_{32}]\mathbb{E}[C_{33}]+\mathbb{E}[z_{31}C_{31}]\mathbb{E}[C_{32}]\mathbb{E}[S_{33}]\} \\ &~+\frac{T_{22}}{T_{32}} \{-\mathbb{E}[S_{31}]\mathbb{E}[z_{32}S_{32}]\mathbb{E}[S_{33}]+\mathbb{E}[S_{31}]\mathbb{E}[z_{32}C_{32}]\mathbb{E}[C_{33}]+\mathbb{E}[C_{31}]\mathbb{E}[z_{32}S_{32}]\mathbb{E}[C_{33}]+\mathbb{E}[C_{31}]\mathbb{E}[z_{32}C_{32}]\mathbb{E}[S_{33}]\} \\ &~+\frac{T_{23}}{T_{33}} \{-\mathbb{E}[S_{31}]\mathbb{E}[S_{32}]\mathbb{E}[z_{33}S_{33}]+\mathbb{E}[S_{31}]\mathbb{E}[C_{32}]\mathbb{E}[z_{33}C_{33}]+\mathbb{E}[C_{31}]\mathbb{E}[S_{32}]\mathbb{E}[z_{33}C_{33}]+\mathbb{E}[C_{31}]\mathbb{E}[C_{32}]\mathbb{E}[z_{33}S_{33}]\} \\ \mathbb{E}&[y^2\cos(\theta)]\\ &=\left(\frac{T_{21}}{T_{31}}\right)^2 \{\mathbb{E}[z_{31}^2C_{31}]\mathbb{E}[C_{32}]\mathbb{E}[C_{33}]-\mathbb{E}[z_{31}^2C_{31}]\mathbb{E}[S_{32}]\mathbb{E}[S_{33}]-\mathbb{E}[z_{31}^2S_{31}]\mathbb{E}[C_{32}]\mathbb{E}[S_{33}]-\mathbb{E}[z_{31}^2S_{31}]\mathbb{E}[S_{32}]\mathbb{E}[C_{33}]\} \\ &~+\left(\frac{T_{22}}{T_{32}}\right)^2 \{\mathbb{E}[C_{31}]\mathbb{E}[z_{32}^2C_{32}]\mathbb{E}[C_{33}]-\mathbb{E}[C_{31}]\mathbb{E}[z_{32}^2S_{32}]\mathbb{E}[S_{33}]-\mathbb{E}[S_{31}]\mathbb{E}[z_{32}^2C_{32}]\mathbb{E}[S_{33}]-\mathbb{E}[S_{31}]\mathbb{E}[z_{32}^2S_{32}]\mathbb{E}[C_{33}]\} \\ &~+\left(\frac{T_{23}}{T_{33}}\right)^2 \{\mathbb{E}[C_{31}]\mathbb{E}[C_{32}]\mathbb{E}[z_{33}^2C_{33}]-\mathbb{E}[C_{31}]\mathbb{E}[S_{32}]\mathbb{E}[z_{33}^2S_{33}]-\mathbb{E}[S_{31}]\mathbb{E}[C_{32}]\mathbb{E}[z_{33}^2S_{33}]-\mathbb{E}[S_{31}]\mathbb{E}[S_{32}]\mathbb{E}[z_{33}^2C_{33}]\} \\ &~+2\frac{T_{21}}{T_{31}}\frac{T_{22}}{T_{32}} \{\mathbb{E}[z_{31}C_{31}]\mathbb{E}[z_{32}C_{32}]\mathbb{E}[C_{33}]-\mathbb{E}[z_{31}C_{31}]\mathbb{E}[z_{32}S_{32}]\mathbb{E}[S_{33}]-\mathbb{E}[z_{31}S_{31}]\mathbb{E}[z_{32}C_{32}]\mathbb{E}[S_{33}]-\mathbb{E}[z_{31}S_{31}]\mathbb{E}[z_{32}S_{32}]\mathbb{E}[C_{33}]\} \\ &~+2\frac{T_{22}}{T_{32}}\frac{T_{23}}{T_{33}} \{\mathbb{E}[C_{31}]\mathbb{E}[z_{32}C_{32}]\mathbb{E}[z_{33}C_{33}]-\mathbb{E}[C_{31}]\mathbb{E}[z_{32}S_{32}]\mathbb{E}[z_{33}S_{33}]-\mathbb{E}[S_{31}]\mathbb{E}[z_{32}C_{32}]\mathbb{E}[z_{33}S_{33}]-\mathbb{E}[S_{31}]\mathbb{E}[z_{32}S_{32}]\mathbb{E}[z_{33}C_{33}]\} \\ &~+2\frac{T_{23}}{T_{33}}\frac{T_{21}}{T_{31}} \{\mathbb{E}[z_{31}C_{31}]\mathbb{E}[C_{32}]\mathbb{E}[z_{33}C_{33}]-\mathbb{E}[z_{31}C_{31}]\mathbb{E}[S_{32}]\mathbb{E}[z_{33}S_{33}]-\mathbb{E}[z_{31}S_{31}]\mathbb{E}[C_{32}]\mathbb{E}[z_{33}S_{33}]-\mathbb{E}[z_{31}S_{31}]\mathbb{E}[S_{32}]\mathbb{E}[z_{33}C_{33}]\} \\ \mathbb{E}&[y^2\sin(\theta)]\\ &=\left(\frac{T_{21}}{T_{31}}\right)^2 \{-\mathbb{E}[z_{31}^2S_{31}]\mathbb{E}[S_{32}]\mathbb{E}[S_{33}]+\mathbb{E}[z_{31}^2S_{31}]\mathbb{E}[C_{32}]\mathbb{E}[C_{33}]+\mathbb{E}[z_{31}^2C_{31}]\mathbb{E}[S_{32}]\mathbb{E}[C_{33}]+\mathbb{E}[z_{31}^2C_{31}]\mathbb{E}[C_{32}]\mathbb{E}[S_{33}]\} \\ &~+\left(\frac{T_{22}}{T_{32}}\right)^2 \{-\mathbb{E}[S_{31}]\mathbb{E}[z_{32}^2S_{32}]\mathbb{E}[S_{33}]+\mathbb{E}[S_{31}]\mathbb{E}[z_{32}^2C_{32}]\mathbb{E}[C_{33}]+\mathbb{E}[C_{31}]\mathbb{E}[z_{32}^2S_{32}]\mathbb{E}[C_{33}]+\mathbb{E}[C_{31}]\mathbb{E}[z_{32}^2C_{32}]\mathbb{E}[S_{33}]\} \\ &~+\left(\frac{T_{23}}{T_{33}}\right)^2 \{-\mathbb{E}[S_{31}]\mathbb{E}[S_{32}]\mathbb{E}[z_{33}^2S_{33}]+\mathbb{E}[S_{31}]\mathbb{E}[C_{32}]\mathbb{E}[z_{33}^2C_{33}]+\mathbb{E}[C_{31}]\mathbb{E}[S_{32}]\mathbb{E}[z_{33}^2C_{33}]+\mathbb{E}[C_{31}]\mathbb{E}[C_{32}]\mathbb{E}[z_{33}^2S_{33}]\} \\ &~+2\frac{T_{21}}{T_{31}}\frac{T_{22}}{T_{32}} \{-\mathbb{E}[z_{31}S_{31}]\mathbb{E}[z_{32}S_{32}]\mathbb{E}[S_{33}]+\mathbb{E}[z_{31}S_{31}]\mathbb{E}[z_{32}C_{32}]\mathbb{E}[C_{33}]+\mathbb{E}[z_{31}C_{31}]\mathbb{E}[z_{32}S_{32}]\mathbb{E}[C_{33}]+\mathbb{E}[z_{31}C_{31}]\mathbb{E}[z_{32}C_{32}]\mathbb{E}[S_{33}]\} \\ &~+2\frac{T_{22}}{T_{32}}\frac{T_{23}}{T_{33}} \{-\mathbb{E}[S_{31}]\mathbb{E}[z_{32}S_{32}]\mathbb{E}[z_{33}S_{33}]+\mathbb{E}[S_{31}]\mathbb{E}[z_{32}C_{32}]\mathbb{E}[z_{33}C_{33}]+\mathbb{E}[C_{31}]\mathbb{E}[z_{32}S_{32}]\mathbb{E}[z_{33}C_{33}]+\mathbb{E}[C_{31}]\mathbb{E}[z_{32}C_{32}]\mathbb{E}[z_{33}S_{33}]\} \\ &~+2\frac{T_{23}}{T_{33}}\frac{T_{21}}{T_{31}} \{-\mathbb{E}[z_{31}S_{31}]\mathbb{E}[S_{32}]\mathbb{E}[z_{33}S_{33}]+\mathbb{E}[z_{31}S_{31}]\mathbb{E}[C_{32}]\mathbb{E}[z_{33}C_{33}]+\mathbb{E}[z_{31}C_{31}]\mathbb{E}[S_{32}]\mathbb{E}[z_{33}C_{33}]+\mathbb{E}[z_{31}C_{31}]\mathbb{E}[C_{32}]\mathbb{E}[z_{33}S_{33}]\} \\ \mathbb{E}&[xy\cos(\theta)]\\ &=\frac{T_{11}}{T_{31}}\frac{T_{21}}{T_{31}} \{\mathbb{E}[z_{31}^2C_{31}]\mathbb{E}[C_{32}]\mathbb{E}[C_{33}]-\mathbb{E}[z_{31}^2C_{31}]\mathbb{E}[S_{32}]\mathbb{E}[S_{33}]-\mathbb{E}[z_{31}^2S_{31}]\mathbb{E}[C_{32}]\mathbb{E}[S_{33}]-\mathbb{E}[z_{31}^2S_{31}]\mathbb{E}[S_{32}]\mathbb{E}[C_{33}]\} \\ &~+\frac{T_{12}}{T_{32}}\frac{T_{22}}{T_{32}} \{\mathbb{E}[C_{31}]\mathbb{E}[z_{32}^2C_{32}]\mathbb{E}[C_{33}]-\mathbb{E}[C_{31}]\mathbb{E}[z_{32}^2S_{32}]\mathbb{E}[S_{33}]-\mathbb{E}[S_{31}]\mathbb{E}[z_{32}^2C_{32}]\mathbb{E}[S_{33}]-\mathbb{E}[S_{31}]\mathbb{E}[z_{32}^2S_{32}]\mathbb{E}[C_{33}]\} \\ &~+\frac{T_{13}}{T_{33}}\frac{T_{23}}{T_{33}} \{\mathbb{E}[C_{31}]\mathbb{E}[C_{32}]\mathbb{E}[z_{33}^2C_{33}]-\mathbb{E}[C_{31}]\mathbb{E}[S_{32}]\mathbb{E}[z_{33}^2S_{33}]-\mathbb{E}[S_{31}]\mathbb{E}[C_{32}]\mathbb{E}[z_{33}^2S_{33}]-\mathbb{E}[S_{31}]\mathbb{E}[S_{32}]\mathbb{E}[z_{33}^2C_{33}]\} \\ &~+\left(\frac{T_{11}}{T_{31}}\frac{T_{22}}{T_{32}}+\frac{T_{12}}{T_{32}}\frac{T_{21}}{T_{31}} \right) \{\mathbb{E}[z_{31}C_{31}]\mathbb{E}[z_{32}C_{32}]\mathbb{E}[C_{33}]-\mathbb{E}[z_{31}C_{31}]\mathbb{E}[z_{32}S_{32}]\mathbb{E}[S_{33}]-\mathbb{E}[z_{31}S_{31}]\mathbb{E}[z_{32}C_{32}]\mathbb{E}[S_{33}]-\mathbb{E}[z_{31}S_{31}]\mathbb{E}[z_{32}S_{32}]\mathbb{E}[C_{33}]\} \\ &~+\left(\frac{T_{12}}{T_{32}}\frac{T_{23}}{T_{33}}+\frac{T_{13}}{T_{33}}\frac{T_{22}}{T_{32}} \right) \{\mathbb{E}[C_{31}]\mathbb{E}[z_{32}C_{32}]\mathbb{E}[z_{33}C_{33}]-\mathbb{E}[C_{31}]\mathbb{E}[z_{32}S_{32}]\mathbb{E}[z_{33}S_{33}]-\mathbb{E}[S_{31}]\mathbb{E}[z_{32}C_{32}]\mathbb{E}[z_{33}S_{33}]-\mathbb{E}[S_{31}]\mathbb{E}[z_{32}S_{32}]\mathbb{E}[z_{33}C_{33}]\} \\ &~+\left(\frac{T_{13}}{T_{31}}\frac{T_{21}}{T_{33}}+\frac{T_{11}}{T_{31}}\frac{T_{23}}{T_{33}} \right) \{\mathbb{E}[z_{31}C_{31}]\mathbb{E}[C_{32}]\mathbb{E}[z_{33}C_{33}]-\mathbb{E}[z_{31}C_{31}]\mathbb{E}[S_{32}]\mathbb{E}[z_{33}S_{33}]-\mathbb{E}[z_{31}S_{31}]\mathbb{E}[C_{32}]\mathbb{E}[z_{33}S_{33}]-\mathbb{E}[z_{31}S_{31}]\mathbb{E}[S_{32}]\mathbb{E}[z_{33}C_{33}]\} \\ \mathbb{E}&[xy\sin(\theta)]\\ &=\frac{T_{11}}{T_{31}}\frac{T_{21}}{T_{31}} \{-\mathbb{E}[z_{31}^2S_{31}]\mathbb{E}[S_{32}]\mathbb{E}[S_{33}]+\mathbb{E}[z_{31}^2S_{31}]\mathbb{E}[C_{32}]\mathbb{E}[C_{33}]+\mathbb{E}[z_{31}^2C_{31}]\mathbb{E}[S_{32}]\mathbb{E}[C_{33}]+\mathbb{E}[z_{31}^2C_{31}]\mathbb{E}[C_{32}]\mathbb{E}[S_{33}]\} \\ &~+\frac{T_{12}}{T_{32}}\frac{T_{22}}{T_{32}} \{-\mathbb{E}[S_{31}]\mathbb{E}[z_{32}^2S_{32}]\mathbb{E}[S_{33}]+\mathbb{E}[S_{31}]\mathbb{E}[z_{32}^2C_{32}]\mathbb{E}[C_{33}]+\mathbb{E}[C_{31}]\mathbb{E}[z_{32}^2S_{32}]\mathbb{E}[C_{33}]+\mathbb{E}[C_{31}]\mathbb{E}[z_{32}^2C_{32}]\mathbb{E}[S_{33}]\} \\ &~+\frac{T_{13}}{T_{33}}\frac{T_{23}}{T_{33}} \{-\mathbb{E}[S_{31}]\mathbb{E}[S_{32}]\mathbb{E}[z_{33}^2S_{33}]+\mathbb{E}[S_{31}]\mathbb{E}[C_{32}]\mathbb{E}[z_{33}^2C_{33}]+\mathbb{E}[C_{31}]\mathbb{E}[S_{32}]\mathbb{E}[z_{33}^2C_{33}]+\mathbb{E}[C_{31}]\mathbb{E}[C_{32}]\mathbb{E}[z_{33}^2S_{33}]\} \\ &~+\left(\frac{T_{11}}{T_{31}}\frac{T_{22}}{T_{32}}+\frac{T_{12}}{T_{32}}\frac{T_{21}}{T_{31}} \right) \{-\mathbb{E}[z_{31}S_{31}]\mathbb{E}[z_{32}S_{32}]\mathbb{E}[S_{33}]+\mathbb{E}[z_{31}S_{31}]\mathbb{E}[z_{32}C_{32}]\mathbb{E}[C_{33}]+\mathbb{E}[z_{31}C_{31}]\mathbb{E}[z_{32}S_{32}]\mathbb{E}[C_{33}]+\mathbb{E}[z_{31}C_{31}]\mathbb{E}[z_{32}C_{32}]\mathbb{E}[S_{33}]\} \\ &~+\left(\frac{T_{12}}{T_{32}}\frac{T_{23}}{T_{33}}+\frac{T_{13}}{T_{33}}\frac{T_{22}}{T_{32}} \right) \{-\mathbb{E}[S_{31}]\mathbb{E}[z_{32}S_{32}]\mathbb{E}[z_{33}S_{33}]+\mathbb{E}[S_{31}]\mathbb{E}[z_{32}C_{32}]\mathbb{E}[z_{33}C_{33}]+\mathbb{E}[C_{31}]\mathbb{E}[z_{32}S_{32}]\mathbb{E}[z_{33}C_{33}]+\mathbb{E}[C_{31}]\mathbb{E}[z_{32}C_{32}]\mathbb{E}[z_{33}S_{33}]\} \\ &~+\left(\frac{T_{13}}{T_{31}}\frac{T_{21}}{T_{33}}+\frac{T_{11}}{T_{31}}\frac{T_{23}}{T_{33}} \right) \{-\mathbb{E}[z_{31}S_{31}]\mathbb{E}[S_{32}]\mathbb{E}[z_{33}S_{33}]+\mathbb{E}[z_{31}S_{31}]\mathbb{E}[C_{32}]\mathbb{E}[z_{33}C_{33}]+\mathbb{E}[z_{31}C_{31}]\mathbb{E}[S_{32}]\mathbb{E}[z_{33}C_{33}]+\mathbb{E}[z_{31}C_{31}]\mathbb{E}[C_{32}]\mathbb{E}[z_{33}S_{33}]\} \\ \mathbb{E}&[x\theta\cos(\theta)]\\ &=\frac{T_{11}}{T_{31}} \{\mathbb{E}[z_{31}^2C_{31}]\mathbb{E}[C_{32}]\mathbb{E}[C_{33}]-\mathbb{E}[z_{31}^2C_{31}]\mathbb{E}[S_{32}]\mathbb{E}[S_{33}]-\mathbb{E}[z_{31}^2S_{31}]\mathbb{E}[C_{32}]\mathbb{E}[S_{33}]-\mathbb{E}[z_{31}^2S_{31}]\mathbb{E}[S_{32}]\mathbb{E}[C_{33}]\} \\ &~+\frac{T_{12}}{T_{32}} \{\mathbb{E}[C_{31}]\mathbb{E}[z_{32}^2C_{32}]\mathbb{E}[C_{33}]-\mathbb{E}[C_{31}]\mathbb{E}[z_{32}^2S_{32}]\mathbb{E}[S_{33}]-\mathbb{E}[S_{31}]\mathbb{E}[z_{32}^2C_{32}]\mathbb{E}[S_{33}]-\mathbb{E}[S_{31}]\mathbb{E}[z_{32}^2S_{32}]\mathbb{E}[C_{33}]\} \\ &~+\frac{T_{13}}{T_{33}} \{\mathbb{E}[C_{31}]\mathbb{E}[C_{32}]\mathbb{E}[z_{33}^2C_{33}]-\mathbb{E}[C_{31}]\mathbb{E}[S_{32}]\mathbb{E}[z_{33}^2S_{33}]-\mathbb{E}[S_{31}]\mathbb{E}[C_{32}]\mathbb{E}[z_{33}^2S_{33}]-\mathbb{E}[S_{31}]\mathbb{E}[S_{32}]\mathbb{E}[z_{33}^2C_{33}]\} \\ &~+\left(\frac{T_{11}}{T_{31}}+\frac{T_{12}}{T_{32}} \right) \{\mathbb{E}[z_{31}C_{31}]\mathbb{E}[z_{32}C_{32}]\mathbb{E}[C_{33}]-\mathbb{E}[z_{31}C_{31}]\mathbb{E}[z_{32}S_{32}]\mathbb{E}[S_{33}]-\mathbb{E}[z_{31}S_{31}]\mathbb{E}[z_{32}C_{32}]\mathbb{E}[S_{33}]-\mathbb{E}[z_{31}S_{31}]\mathbb{E}[z_{32}S_{32}]\mathbb{E}[C_{33}]\} \\ &~+\left(\frac{T_{12}}{T_{32}}+\frac{T_{13}}{T_{33}} \right) \{\mathbb{E}[C_{31}]\mathbb{E}[z_{32}C_{32}]\mathbb{E}[z_{33}C_{33}]-\mathbb{E}[C_{31}]\mathbb{E}[z_{32}S_{32}]\mathbb{E}[z_{33}S_{33}]-\mathbb{E}[S_{31}]\mathbb{E}[z_{32}C_{32}]\mathbb{E}[z_{33}S_{33}]-\mathbb{E}[S_{31}]\mathbb{E}[z_{32}S_{32}]\mathbb{E}[z_{33}C_{33}]\} \\ &~+\left(\frac{T_{13}}{T_{31}}+\frac{T_{11}}{T_{31}} \right) \{\mathbb{E}[z_{31}C_{31}]\mathbb{E}[C_{32}]\mathbb{E}[z_{33}C_{33}]-\mathbb{E}[z_{31}C_{31}]\mathbb{E}[S_{32}]\mathbb{E}[z_{33}S_{33}]-\mathbb{E}[z_{31}S_{31}]\mathbb{E}[C_{32}]\mathbb{E}[z_{33}S_{33}]-\mathbb{E}[z_{31}S_{31}]\mathbb{E}[S_{32}]\mathbb{E}[z_{33}C_{33}]\} \\ \mathbb{E}&[x\theta\sin(\theta)]\\ &=\frac{T_{11}}{T_{31}} \{-\mathbb{E}[z_{31}^2S_{31}]\mathbb{E}[S_{32}]\mathbb{E}[S_{33}]+\mathbb{E}[z_{31}^2S_{31}]\mathbb{E}[C_{32}]\mathbb{E}[C_{33}]+\mathbb{E}[z_{31}^2C_{31}]\mathbb{E}[S_{32}]\mathbb{E}[C_{33}]+\mathbb{E}[z_{31}^2C_{31}]\mathbb{E}[C_{32}]\mathbb{E}[S_{33}]\} \\ &~+\frac{T_{12}}{T_{32}} \{-\mathbb{E}[S_{31}]\mathbb{E}[z_{32}^2S_{32}]\mathbb{E}[S_{33}]+\mathbb{E}[S_{31}]\mathbb{E}[z_{32}^2C_{32}]\mathbb{E}[C_{33}]+\mathbb{E}[C_{31}]\mathbb{E}[z_{32}^2S_{32}]\mathbb{E}[C_{33}]+\mathbb{E}[C_{31}]\mathbb{E}[z_{32}^2C_{32}]\mathbb{E}[S_{33}]\} \\ &~+\frac{T_{13}}{T_{33}} \{-\mathbb{E}[S_{31}]\mathbb{E}[S_{32}]\mathbb{E}[z_{33}^2S_{33}]+\mathbb{E}[S_{31}]\mathbb{E}[C_{32}]\mathbb{E}[z_{33}^2C_{33}]+\mathbb{E}[C_{31}]\mathbb{E}[S_{32}]\mathbb{E}[z_{33}^2C_{33}]+\mathbb{E}[C_{31}]\mathbb{E}[C_{32}]\mathbb{E}[z_{33}^2S_{33}]\} \\ &~+\left(\frac{T_{11}}{T_{31}}+\frac{T_{12}}{T_{32}} \right) \{-\mathbb{E}[z_{31}S_{31}]\mathbb{E}[z_{32}S_{32}]\mathbb{E}[S_{33}]+\mathbb{E}[z_{31}S_{31}]\mathbb{E}[z_{32}C_{32}]\mathbb{E}[C_{33}]+\mathbb{E}[z_{31}C_{31}]\mathbb{E}[z_{32}S_{32}]\mathbb{E}[C_{33}]+\mathbb{E}[z_{31}C_{31}]\mathbb{E}[z_{32}C_{32}]\mathbb{E}[S_{33}]\} \\ &~+\left(\frac{T_{12}}{T_{32}}+\frac{T_{13}}{T_{33}} \right) \{-\mathbb{E}[S_{31}]\mathbb{E}[z_{32}S_{32}]\mathbb{E}[z_{33}S_{33}]+\mathbb{E}[S_{31}]\mathbb{E}[z_{32}C_{32}]\mathbb{E}[z_{33}C_{33}]+\mathbb{E}[C_{31}]\mathbb{E}[z_{32}S_{32}]\mathbb{E}[z_{33}C_{33}]+\mathbb{E}[C_{31}]\mathbb{E}[z_{32}C_{32}]\mathbb{E}[z_{33}S_{33}]\} \\ &~+\left(\frac{T_{13}}{T_{31}}+\frac{T_{11}}{T_{31}} \right) \{-\mathbb{E}[z_{31}S_{31}]\mathbb{E}[S_{32}]\mathbb{E}[z_{33}S_{33}]+\mathbb{E}[z_{31}S_{31}]\mathbb{E}[C_{32}]\mathbb{E}[z_{33}C_{33}]+\mathbb{E}[z_{31}C_{31}]\mathbb{E}[S_{32}]\mathbb{E}[z_{33}C_{33}]+\mathbb{E}[z_{31}C_{31}]\mathbb{E}[C_{32}]\mathbb{E}[z_{33}S_{33}]\} \\ \mathbb{E}&[y\theta\cos(\theta)]\\ &=\frac{T_{21}}{T_{31}} \{\mathbb{E}[z_{31}^2C_{31}]\mathbb{E}[C_{32}]\mathbb{E}[C_{33}]-\mathbb{E}[z_{31}^2C_{31}]\mathbb{E}[S_{32}]\mathbb{E}[S_{33}]-\mathbb{E}[z_{31}^2S_{31}]\mathbb{E}[C_{32}]\mathbb{E}[S_{33}]-\mathbb{E}[z_{31}^2S_{31}]\mathbb{E}[S_{32}]\mathbb{E}[C_{33}]\} \\ &~+\frac{T_{22}}{T_{32}} \{\mathbb{E}[C_{31}]\mathbb{E}[z_{32}^2C_{32}]\mathbb{E}[C_{33}]-\mathbb{E}[C_{31}]\mathbb{E}[z_{32}^2S_{32}]\mathbb{E}[S_{33}]-\mathbb{E}[S_{31}]\mathbb{E}[z_{32}^2C_{32}]\mathbb{E}[S_{33}]-\mathbb{E}[S_{31}]\mathbb{E}[z_{32}^2S_{32}]\mathbb{E}[C_{33}]\} \\ &~+\frac{T_{23}}{T_{33}} \{\mathbb{E}[C_{31}]\mathbb{E}[C_{32}]\mathbb{E}[z_{33}^2C_{33}]-\mathbb{E}[C_{31}]\mathbb{E}[S_{32}]\mathbb{E}[z_{33}^2S_{33}]-\mathbb{E}[S_{31}]\mathbb{E}[C_{32}]\mathbb{E}[z_{33}^2S_{33}]-\mathbb{E}[S_{31}]\mathbb{E}[S_{32}]\mathbb{E}[z_{33}^2C_{33}]\} \\ &~+\left(\frac{T_{21}}{T_{31}}+\frac{T_{22}}{T_{32}} \right) \{\mathbb{E}[z_{31}C_{31}]\mathbb{E}[z_{32}C_{32}]\mathbb{E}[C_{33}]-\mathbb{E}[z_{31}C_{31}]\mathbb{E}[z_{32}S_{32}]\mathbb{E}[S_{33}]-\mathbb{E}[z_{31}S_{31}]\mathbb{E}[z_{32}C_{32}]\mathbb{E}[S_{33}]-\mathbb{E}[z_{31}S_{31}]\mathbb{E}[z_{32}S_{32}]\mathbb{E}[C_{33}]\} \\ &~+\left(\frac{T_{22}}{T_{32}}+\frac{T_{23}}{T_{33}} \right) \{\mathbb{E}[C_{31}]\mathbb{E}[z_{32}C_{32}]\mathbb{E}[z_{33}C_{33}]-\mathbb{E}[C_{31}]\mathbb{E}[z_{32}S_{32}]\mathbb{E}[z_{33}S_{33}]-\mathbb{E}[S_{31}]\mathbb{E}[z_{32}C_{32}]\mathbb{E}[z_{33}S_{33}]-\mathbb{E}[S_{31}]\mathbb{E}[z_{32}S_{32}]\mathbb{E}[z_{33}C_{33}]\} \\ &~+\left(\frac{T_{23}}{T_{31}}+\frac{T_{21}}{T_{31}} \right) \{\mathbb{E}[z_{31}C_{31}]\mathbb{E}[C_{32}]\mathbb{E}[z_{33}C_{33}]-\mathbb{E}[z_{31}C_{31}]\mathbb{E}[S_{32}]\mathbb{E}[z_{33}S_{33}]-\mathbb{E}[z_{31}S_{31}]\mathbb{E}[C_{32}]\mathbb{E}[z_{33}S_{33}]-\mathbb{E}[z_{31}S_{31}]\mathbb{E}[S_{32}]\mathbb{E}[z_{33}C_{33}]\} \\ \mathbb{E}&[y\theta\sin(\theta)]\\ &=\frac{T_{21}}{T_{31}} \{-\mathbb{E}[z_{31}^2S_{31}]\mathbb{E}[S_{32}]\mathbb{E}[S_{33}]+\mathbb{E}[z_{31}^2S_{31}]\mathbb{E}[C_{32}]\mathbb{E}[C_{33}]+\mathbb{E}[z_{31}^2C_{31}]\mathbb{E}[S_{32}]\mathbb{E}[C_{33}]+\mathbb{E}[z_{31}^2C_{31}]\mathbb{E}[C_{32}]\mathbb{E}[S_{33}]\} \\ &~+\frac{T_{22}}{T_{32}} \{-\mathbb{E}[S_{31}]\mathbb{E}[z_{32}^2S_{32}]\mathbb{E}[S_{33}]+\mathbb{E}[S_{31}]\mathbb{E}[z_{32}^2C_{32}]\mathbb{E}[C_{33}]+\mathbb{E}[C_{31}]\mathbb{E}[z_{32}^2S_{32}]\mathbb{E}[C_{33}]+\mathbb{E}[C_{31}]\mathbb{E}[z_{32}^2C_{32}]\mathbb{E}[S_{33}]\} \\ &~+\frac{T_{23}}{T_{33}} \{-\mathbb{E}[S_{31}]\mathbb{E}[S_{32}]\mathbb{E}[z_{33}^2S_{33}]+\mathbb{E}[S_{31}]\mathbb{E}[C_{32}]\mathbb{E}[z_{33}^2C_{33}]+\mathbb{E}[C_{31}]\mathbb{E}[S_{32}]\mathbb{E}[z_{33}^2C_{33}]+\mathbb{E}[C_{31}]\mathbb{E}[C_{32}]\mathbb{E}[z_{33}^2S_{33}]\} \\ &~+\left(\frac{T_{21}}{T_{31}}+\frac{T_{22}}{T_{32}} \right) \{-\mathbb{E}[z_{31}S_{31}]\mathbb{E}[z_{32}S_{32}]\mathbb{E}[S_{33}]+\mathbb{E}[z_{31}S_{31}]\mathbb{E}[z_{32}C_{32}]\mathbb{E}[C_{33}]+\mathbb{E}[z_{31}C_{31}]\mathbb{E}[z_{32}S_{32}]\mathbb{E}[C_{33}]+\mathbb{E}[z_{31}C_{31}]\mathbb{E}[z_{32}C_{32}]\mathbb{E}[S_{33}]\} \\ &~+\left(\frac{T_{22}}{T_{32}}+\frac{T_{23}}{T_{33}} \right) \{-\mathbb{E}[S_{31}]\mathbb{E}[z_{32}S_{32}]\mathbb{E}[z_{33}S_{33}]+\mathbb{E}[S_{31}]\mathbb{E}[z_{32}C_{32}]\mathbb{E}[z_{33}C_{33}]+\mathbb{E}[C_{31}]\mathbb{E}[z_{32}S_{32}]\mathbb{E}[z_{33}C_{33}]+\mathbb{E}[C_{31}]\mathbb{E}[z_{32}C_{32}]\mathbb{E}[z_{33}S_{33}]\} \\ &~+\left(\frac{T_{23}}{T_{31}}+\frac{T_{21}}{T_{31}} \right) \{-\mathbb{E}[z_{31}S_{31}]\mathbb{E}[S_{32}]\mathbb{E}[z_{33}S_{33}]+\mathbb{E}[z_{31}S_{31}]\mathbb{E}[C_{32}]\mathbb{E}[z_{33}C_{33}]+\mathbb{E}[z_{31}C_{31}]\mathbb{E}[S_{32}]\mathbb{E}[z_{33}C_{33}]+\mathbb{E}[z_{31}C_{31}]\mathbb{E}[C_{32}]\mathbb{E}[z_{33}S_{33}]\} \\ \end{align*}

計算(MATLAB)とモンテカルロ法での確認

2変数の場合

確率変数を下記のように設定する．

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

結果は下記で，正しく計算できているといえる．

\mathbb{E}[x\theta] \mathbb{E}[x\cos(\theta)] \mathbb{E}[x\sin(\theta)] \mathbb{E}[x\cos(\theta)\sin(\theta)]
解析解 11.97 2.85 7.24 1.26
モンテカルロ 11.97 2.85 7.24 1.26
差( \times 10^{-4} ) -7.19 4.36 -0.19 3.77
コード(厳密解)
% 確率変数の定義
mu_x = [10.0;pi/3];
Sigma_x = [5.0,1.5;1.5,pi/6];
% 直交行列と固有値の計算
[T,lambda] = eig(Sigma_x);
T11 = T(1,1);
T12 = T(1,2);
T21 = T(2,1);
T22 = T(2,2);
lambda1 = lambda(1,1);
lambda2 = lambda(2,2);
% mu_yの計算
mu_y = T'*mu_x;
mu_y1 = mu_y(1);
mu_y2 = mu_y(2);

% 確率分布のクラス定義
y1 = GaussianMTPM(mu_y1,lambda1);
y2 = GaussianMTPM(mu_y2,lambda2);
l1 = GaussianMTPM(T21*mu_y1,T21^2*lambda1);
l2 = GaussianMTPM(T22*mu_y2,T22^2*lambda2);
m1 = GaussianMTPM(2*T21*mu_y1,(2*T21)^2*lambda1);
m2 = GaussianMTPM(2*T22*mu_y2,(2*T22)^2*(lambda2));

% 計算
% E[x*theta]
E_xth = T11*T21*y1.X2+(T11*T22+T12*T21)*y1.X*y2.X+T12*T22*y2.X2
% E[x*cos(theta)]
E_xC = T11/T21*(l1.XCosX*l2.CosX-l1.XSinX*l2.SinX)+T12/T22*(l1.CosX*l2.XCosX-l1.SinX*l2.XSinX)
% E[x*sin(theta)]
E_xS = T11/T21*(l1.XSinX*l2.CosX+l1.XCosX*l2.SinX)+T12/T22*(l1.SinX*l2.XCosX+l1.CosX*l2.XSinX)
% E[x*cos(theta)*sin(theta)]
E_xCS = T11/(4*T21)*(m1.XSinX*m2.CosX+m1.XCosX*m2.SinX)+T12/(4*T22)*(m1.SinX*m2.XCosX+m1.CosX*m2.XSinX)
コード（モンテカルロ法）
rng('default');
rng(1);
x = mvnrnd(mu_x,Sigma_x,1e8);
M_xth=mean(x(:,1).*x(:,2))
M_xC=mean(x(:,1).*cos(x(:,2)))
M_xS=mean(x(:,1).*sin(x(:,2)))
M_xCS=mean(x(:,1).*cos(x(:,2)).*sin(x(:,2)))
コード（結果の表示用）

e1 = E_xth-M_xth
e2 = E_xC-M_xC
e3 = E_xS-M_xS
e4 = E_xCS-M_xCS

s1 = sprintf('| | \\\\mathbb{E}[x\\\\theta] | \\\\mathbb{E}[x\\\\cos(\\\\theta)] | \\\\mathbb{E}[x\\\\sin(\\\\theta)] | \\\\mathbb{E}[x\\\\cos(\\\\theta)\\\\sin(\\\\theta)] |\n');
s2 = sprintf('| ---- | ---- | ---- | ---- | ---- |\n');
s3 = sprintf('|解析解|%.2f|%.2f|%.2f|%.2f|\n',E_xth,E_xC,E_xS,E_xCS);
s4 = sprintf('|モンテカルロ|%.2f|%.2f|%.2f|%.2f|\n',M_xth,M_xC,M_xS,M_xCS);
s5 = sprintf('|差( \\\\times 10^{-4} )|%.2f|%.2f|%.2f|%.2f|\n',e11e4,e21e4,e31e4,e41e4);
fprintf([s1,s2,s3,s4,s5])

3変数の場合

確率変数を下記のように設定する．

\begin{align*} & \begin{pmatrix} x\\ \theta \end{pmatrix} \sim \mathcal{N} \left( \begin{pmatrix} {10.0}\\5.0\\{\pi/3} \end{pmatrix}, \begin{pmatrix} 5.0&0.5&0.5\\ 0.5&2.0&0.3\\0.5&0.3&\pi/10 \end{pmatrix} \right) \end{align*}

結果は下記で，こちらも正しく計算できているといえる．

\mathbb{E}[x\cos(\theta)] \mathbb{E}[x\sin(\theta)] \mathbb{E}[y\cos(\theta)] \mathbb{E}[y\sin(\theta)] \mathbb{E}[xy\cos(\theta)] \mathbb{E}[xy\sin(\theta)] \mathbb{E}[x^2\cos(\theta)] \mathbb{E}[x^2\sin(\theta)] \mathbb{E}[y^2\cos(\theta)] \mathbb{E}[y^2\sin(\theta)] \mathbb{E}[x\theta\cos(\theta)] \mathbb{E}[x\theta\sin(\theta)] \mathbb{E}[y\theta\cos(\theta)] \mathbb{E}[y\theta\sin(\theta)]
解析解 3.90 7.62 1.91 3.83 17.44 39.62 36.51 80.32 9.28 21.20 1.91 9.57 0.93 4.83
モンテカルロ 3.90 7.62 1.91 3.83 17.44 39.62 36.51 80.32 9.28 21.20 1.91 9.57 0.93 4.83
差( \times 10^{-4} ) -1.15 -2.25 0.40 -1.06 -2.34 -13.32 -23.62 -22.44 4.05 -10.03 -2.89 -1.37 -0.53 -1.49
コード(厳密解)
% 確率変数の定義
mu_x = [10.0;5.0;pi/3]
Sigma_x = [3,0.5,0.5;0.5,2.0,0.3;0.5,0.3,pi/10;]

% 直交行列と固有値の計算
[T,lambda] = eig(Sigma_x);
T11 = T(1,1);
T12 = T(1,2);
T13 = T(1,3);
T21 = T(2,1);
T22 = T(2,2);
T23 = T(2,3);
T31 = T(3,1);
T32 = T(3,2);
T33 = T(3,3);

lambda1 = lambda(1,1);
lambda2 = lambda(2,2);  
lambda3 = lambda(3,3);

mu_y = T'*mu_x;
mu_y1 = mu_y(1);
mu_y2 = mu_y(2);
mu_y3 = mu_y(3);

% 確率分布のクラス定義
y1 = GaussianMTPM(mu_y1,lambda1);
y2 = GaussianMTPM(mu_y2,lambda2);
y3 = GaussianMTPM(mu_y3,lambda3);
z31 = GaussianMTPM(T31*mu_y1,(T31^2)*lambda1);
z32 = GaussianMTPM(T32*mu_y2,(T32^2)*lambda2);
z33 = GaussianMTPM(T33*mu_y3,(T33^2)*lambda3);

% calc E[(x*cos(theta))]
E_xC =  (T11/T31)*(z31.XCosX*z32.CosX*z33.CosX-z31.XCosX*z32.SinX*z33.SinX-z31.XSinX*z32.CosX*z33.SinX-z31.XSinX*z32.SinX*z33.CosX)...
    +(T12/T32)*(z31.CosX*z32.XCosX*z33.CosX-z31.CosX*z32.XSinX*z33.SinX-z31.SinX*z32.XCosX*z33.SinX-z31.SinX*z32.XSinX*z33.CosX)...
    +(T13/T33)*(z31.CosX*z32.CosX*z33.XCosX-z31.CosX*z32.SinX*z33.XSinX-z31.SinX*z32.CosX*z33.XSinX-z31.SinX*z32.SinX*z33.XCosX)

% calc E[(x*sin(theta))]
E_xS =  (T11/T31)*(-z31.XSinX*z32.SinX*z33.SinX+z31.XSinX*z32.CosX*z33.CosX+z31.XCosX*z32.SinX*z33.CosX+z31.XCosX*z32.CosX*z33.SinX)...
    +(T12/T32)*(-z31.SinX*z32.XSinX*z33.SinX+z31.SinX*z32.XCosX*z33.CosX+z31.CosX*z32.XSinX*z33.CosX+z31.CosX*z32.XCosX*z33.SinX)...
    +(T13/T33)*(-z31.SinX*z32.SinX*z33.XSinX+z31.SinX*z32.CosX*z33.XCosX+z31.CosX*z32.SinX*z33.XCosX+z31.CosX*z32.CosX*z33.XSinX)

% calc E[(y*cos(theta))]
E_yC =  (T21/T31)*(z31.XCosX*z32.CosX*z33.CosX-z31.XCosX*z32.SinX*z33.SinX-z31.XSinX*z32.CosX*z33.SinX-z31.XSinX*z32.SinX*z33.CosX)...
    +(T22/T32)*(z31.CosX*z32.XCosX*z33.CosX-z31.CosX*z32.XSinX*z33.SinX-z31.SinX*z32.XCosX*z33.SinX-z31.SinX*z32.XSinX*z33.CosX)...
    +(T23/T33)*(z31.CosX*z32.CosX*z33.XCosX-z31.CosX*z32.SinX*z33.XSinX-z31.SinX*z32.CosX*z33.XSinX-z31.SinX*z32.SinX*z33.XCosX)

% calc E[(y*sin(theta))]
% sin:: (-z31.SinX*z32.SinX*z33.SinX+z31.SinX*z32.CosX*z33.CosX+z31.CosX*z32.SinX*z33.CosX+z31.CosX*z32.CosX*z33.SinX)
E_yS =  (T21/T31)*(-z31.XSinX*z32.SinX*z33.SinX+z31.XSinX*z32.CosX*z33.CosX+z31.XCosX*z32.SinX*z33.CosX+z31.XCosX*z32.CosX*z33.SinX)...
    +(T22/T32)*(-z31.SinX*z32.XSinX*z33.SinX+z31.SinX*z32.XCosX*z33.CosX+z31.CosX*z32.XSinX*z33.CosX+z31.CosX*z32.XCosX*z33.SinX)...
    +(T23/T33)*(-z31.SinX*z32.SinX*z33.XSinX+z31.SinX*z32.CosX*z33.XCosX+z31.CosX*z32.SinX*z33.XCosX+z31.CosX*z32.CosX*z33.XSinX)

% calc E[(x*x*cos(theta))]
E_x2C =  (T11/T31)^2*(z31.X2CosX*z32.CosX*z33.CosX-z31.X2CosX*z32.SinX*z33.SinX-z31.X2SinX*z32.CosX*z33.SinX-z31.X2SinX*z32.SinX*z33.CosX)...
    +(T12/T32)^2*(z31.CosX*z32.X2CosX*z33.CosX-z31.CosX*z32.X2SinX*z33.SinX-z31.SinX*z32.X2CosX*z33.SinX-z31.SinX*z32.X2SinX*z33.CosX)...
    +(T13/T33)^2*(z31.CosX*z32.CosX*z33.X2CosX-z31.CosX*z32.SinX*z33.X2SinX-z31.SinX*z32.CosX*z33.X2SinX-z31.SinX*z32.SinX*z33.X2CosX)...
    +2*(T11/T31)*(T12/T32)*(z31.XCosX*z32.XCosX*z33.CosX-z31.XCosX*z32.XSinX*z33.SinX-z31.XSinX*z32.XCosX*z33.SinX-z31.XSinX*z32.XSinX*z33.CosX)...
    +2*(T12/T32)*(T13/T33)*(z31.CosX*z32.XCosX*z33.XCosX-z31.CosX*z32.XSinX*z33.XSinX-z31.SinX*z32.XCosX*z33.XSinX-z31.SinX*z32.XSinX*z33.XCosX)...
    +2*(T13/T33)*(T11/T31)*(z31.XCosX*z32.CosX*z33.XCosX-z31.XCosX*z32.SinX*z33.XSinX-z31.XSinX*z32.CosX*z33.XSinX-z31.XSinX*z32.SinX*z33.XCosX)

% calc E[(x*x*sin(theta))]
% sin:: (-z31.SinX*z32.SinX*z33.SinX+z31.SinX*z32.CosX*z33.CosX+z31.CosX*z32.SinX*z33.CosX+z31.CosX*z32.CosX*z33.SinX)
E_x2S =  (T11/T31)^2*(-z31.X2SinX*z32.SinX*z33.SinX+z31.X2SinX*z32.CosX*z33.CosX+z31.X2CosX*z32.SinX*z33.CosX+z31.X2CosX*z32.CosX*z33.SinX)...
    +(T12/T32)^2*(-z31.SinX*z32.X2SinX*z33.SinX+z31.SinX*z32.X2CosX*z33.CosX+z31.CosX*z32.X2SinX*z33.CosX+z31.CosX*z32.X2CosX*z33.SinX)...
    +(T13/T33)^2*(-z31.SinX*z32.SinX*z33.X2SinX+z31.SinX*z32.CosX*z33.X2CosX+z31.CosX*z32.SinX*z33.X2CosX+z31.CosX*z32.CosX*z33.X2SinX)...
    +2*(T11/T31)*(T12/T32)*(-z31.XSinX*z32.XSinX*z33.SinX+z31.XSinX*z32.XCosX*z33.CosX+z31.XCosX*z32.XSinX*z33.CosX+z31.XCosX*z32.XCosX*z33.SinX)...
    +2*(T12/T32)*(T13/T33)*(-z31.SinX*z32.XSinX*z33.XSinX+z31.SinX*z32.XCosX*z33.XCosX+z31.CosX*z32.XSinX*z33.XCosX+z31.CosX*z32.XCosX*z33.XSinX)...
    +2*(T13/T33)*(T11/T31)*(-z31.XSinX*z32.SinX*z33.XSinX+z31.XSinX*z32.CosX*z33.XCosX+z31.XCosX*z32.SinX*z33.XCosX+z31.XCosX*z32.CosX*z33.XSinX)

% calc E[(x*y*cos(theta))]
E_xyC =  (T11/T31)*(T21/T31)*(z31.X2CosX*z32.CosX*z33.CosX-z31.X2CosX*z32.SinX*z33.SinX-z31.X2SinX*z32.CosX*z33.SinX-z31.X2SinX*z32.SinX*z33.CosX)...
    +(T12/T32)*(T22/T32)*(z31.CosX*z32.X2CosX*z33.CosX-z31.CosX*z32.X2SinX*z33.SinX-z31.SinX*z32.X2CosX*z33.SinX-z31.SinX*z32.X2SinX*z33.CosX)...
    +(T13/T33)*(T23/T33)*(z31.CosX*z32.CosX*z33.X2CosX-z31.CosX*z32.SinX*z33.X2SinX-z31.SinX*z32.CosX*z33.X2SinX-z31.SinX*z32.SinX*z33.X2CosX)...
    +((T11/T31)*(T22/T32)+(T12/T32)*(T21/T31))*(z31.XCosX*z32.XCosX*z33.CosX-z31.XCosX*z32.XSinX*z33.SinX-z31.XSinX*z32.XCosX*z33.SinX-z31.XSinX*z32.XSinX*z33.CosX)...
    +((T12/T32)*(T23/T33)+(T13/T33)*(T22/T32))*(z31.CosX*z32.XCosX*z33.XCosX-z31.CosX*z32.XSinX*z33.XSinX-z31.SinX*z32.XCosX*z33.XSinX-z31.SinX*z32.XSinX*z33.XCosX)...
    +((T13/T33)*(T21/T31)+(T11/T31)*(T23/T33))*(z31.XCosX*z32.CosX*z33.XCosX-z31.XCosX*z32.SinX*z33.XSinX-z31.XSinX*z32.CosX*z33.XSinX-z31.XSinX*z32.SinX*z33.XCosX)

% calc E[(x*y*sin(theta))]
E_xyS =  (T11/T31)*(T21/T31)*(-z31.X2SinX*z32.SinX*z33.SinX+z31.X2SinX*z32.CosX*z33.CosX+z31.X2CosX*z32.SinX*z33.CosX+z31.X2CosX*z32.CosX*z33.SinX)...
    +(T12/T32)*(T22/T32)*(-z31.SinX*z32.X2SinX*z33.SinX+z31.SinX*z32.X2CosX*z33.CosX+z31.CosX*z32.X2SinX*z33.CosX+z31.CosX*z32.X2CosX*z33.SinX)...
    +(T13/T33)*(T23/T33)*(-z31.SinX*z32.SinX*z33.X2SinX+z31.SinX*z32.CosX*z33.X2CosX+z31.CosX*z32.SinX*z33.X2CosX+z31.CosX*z32.CosX*z33.X2SinX)...
    +((T11/T31)*(T22/T32)+(T12/T32)*(T21/T31))*(-z31.XSinX*z32.XSinX*z33.SinX+z31.XSinX*z32.XCosX*z33.CosX+z31.XCosX*z32.XSinX*z33.CosX+z31.XCosX*z32.XCosX*z33.SinX)...
    +((T12/T32)*(T23/T33)+(T13/T33)*(T22/T32))*(-z31.SinX*z32.XSinX*z33.XSinX+z31.SinX*z32.XCosX*z33.XCosX+z31.CosX*z32.XSinX*z33.XCosX+z31.CosX*z32.XCosX*z33.XSinX)...
    +((T13/T33)*(T21/T31)+(T11/T31)*(T23/T33))*(-z31.XSinX*z32.SinX*z33.XSinX+z31.XSinX*z32.CosX*z33.XCosX+z31.XCosX*z32.SinX*z33.XCosX+z31.XCosX*z32.CosX*z33.XSinX)

% calc E[(y*y*cos(theta))]
E_y2C =  (T21/T31)^2*(z31.X2CosX*z32.CosX*z33.CosX-z31.X2CosX*z32.SinX*z33.SinX-z31.X2SinX*z32.CosX*z33.SinX-z31.X2SinX*z32.SinX*z33.CosX)...
    +(T22/T32)^2*(z31.CosX*z32.X2CosX*z33.CosX-z31.CosX*z32.X2SinX*z33.SinX-z31.SinX*z32.X2CosX*z33.SinX-z31.SinX*z32.X2SinX*z33.CosX)...
    +(T23/T33)^2*(z31.CosX*z32.CosX*z33.X2CosX-z31.CosX*z32.SinX*z33.X2SinX-z31.SinX*z32.CosX*z33.X2SinX-z31.SinX*z32.SinX*z33.X2CosX)...
    +2*(T21/T31)*(T22/T32)*(z31.XCosX*z32.XCosX*z33.CosX-z31.XCosX*z32.XSinX*z33.SinX-z31.XSinX*z32.XCosX*z33.SinX-z31.XSinX*z32.XSinX*z33.CosX)...
    +2*(T22/T32)*(T23/T33)*(z31.CosX*z32.XCosX*z33.XCosX-z31.CosX*z32.XSinX*z33.XSinX-z31.SinX*z32.XCosX*z33.XSinX-z31.SinX*z32.XSinX*z33.XCosX)...
    +2*(T23/T33)*(T21/T31)*(z31.XCosX*z32.CosX*z33.XCosX-z31.XCosX*z32.SinX*z33.XSinX-z31.XSinX*z32.CosX*z33.XSinX-z31.XSinX*z32.SinX*z33.XCosX)

% calc E[(y*y*sin(theta))]
E_y2S =  (T21/T31)^2*(-z31.X2SinX*z32.SinX*z33.SinX+z31.X2SinX*z32.CosX*z33.CosX+z31.X2CosX*z32.SinX*z33.CosX+z31.X2CosX*z32.CosX*z33.SinX)...
    +(T22/T32)^2*(-z31.SinX*z32.X2SinX*z33.SinX+z31.SinX*z32.X2CosX*z33.CosX+z31.CosX*z32.X2SinX*z33.CosX+z31.CosX*z32.X2CosX*z33.SinX)...
    +(T23/T33)^2*(-z31.SinX*z32.SinX*z33.X2SinX+z31.SinX*z32.CosX*z33.X2CosX+z31.CosX*z32.SinX*z33.X2CosX+z31.CosX*z32.CosX*z33.X2SinX)...
    +2*(T21/T31)*(T22/T32)*(-z31.XSinX*z32.XSinX*z33.SinX+z31.XSinX*z32.XCosX*z33.CosX+z31.XCosX*z32.XSinX*z33.CosX+z31.XCosX*z32.XCosX*z33.SinX)...
    +2*(T22/T32)*(T23/T33)*(-z31.SinX*z32.XSinX*z33.XSinX+z31.SinX*z32.XCosX*z33.XCosX+z31.CosX*z32.XSinX*z33.XCosX+z31.CosX*z32.XCosX*z33.XSinX)...
    +2*(T23/T33)*(T21/T31)*(-z31.XSinX*z32.SinX*z33.XSinX+z31.XSinX*z32.CosX*z33.XCosX+z31.XCosX*z32.SinX*z33.XCosX+z31.XCosX*z32.CosX*z33.XSinX)

% calc E[(x*theta*cos(theta))]
E_xthC =  (T11/T31)*(z31.X2CosX*z32.CosX*z33.CosX-z31.X2CosX*z32.SinX*z33.SinX-z31.X2SinX*z32.CosX*z33.SinX-z31.X2SinX*z32.SinX*z33.CosX)...
    +(T12/T32)*(z31.CosX*z32.X2CosX*z33.CosX-z31.CosX*z32.X2SinX*z33.SinX-z31.SinX*z32.X2CosX*z33.SinX-z31.SinX*z32.X2SinX*z33.CosX)...
    +(T13/T33)*(z31.CosX*z32.CosX*z33.X2CosX-z31.CosX*z32.SinX*z33.X2SinX-z31.SinX*z32.CosX*z33.X2SinX-z31.SinX*z32.SinX*z33.X2CosX)...
    +((T11/T31)+(T12/T32))*(z31.XCosX*z32.XCosX*z33.CosX-z31.XCosX*z32.XSinX*z33.SinX-z31.XSinX*z32.XCosX*z33.SinX-z31.XSinX*z32.XSinX*z33.CosX)...
    +((T12/T32)+(T13/T33))*(z31.CosX*z32.XCosX*z33.XCosX-z31.CosX*z32.XSinX*z33.XSinX-z31.SinX*z32.XCosX*z33.XSinX-z31.SinX*z32.XSinX*z33.XCosX)...
    +((T11/T31)+(T13/T33))*(z31.XCosX*z32.CosX*z33.XCosX-z31.XCosX*z32.SinX*z33.XSinX-z31.XSinX*z32.CosX*z33.XSinX-z31.XSinX*z32.SinX*z33.XCosX)

% calc E[(x*theta*sin(theta))]
E_xthS =  (T11/T31)*(-z31.X2SinX*z32.SinX*z33.SinX+z31.X2SinX*z32.CosX*z33.CosX+z31.X2CosX*z32.SinX*z33.CosX+z31.X2CosX*z32.CosX*z33.SinX)...
    +(T12/T32)*(-z31.SinX*z32.X2SinX*z33.SinX+z31.SinX*z32.X2CosX*z33.CosX+z31.CosX*z32.X2SinX*z33.CosX+z31.CosX*z32.X2CosX*z33.SinX)...
    +(T13/T33)*(-z31.SinX*z32.SinX*z33.X2SinX+z31.SinX*z32.CosX*z33.X2CosX+z31.CosX*z32.SinX*z33.X2CosX+z31.CosX*z32.CosX*z33.X2SinX)...
    +((T11/T31)+(T12/T32))*(-z31.XSinX*z32.XSinX*z33.SinX+z31.XSinX*z32.XCosX*z33.CosX+z31.XCosX*z32.XSinX*z33.CosX+z31.XCosX*z32.XCosX*z33.SinX)...
    +((T12/T32)+(T13/T33))*(-z31.SinX*z32.XSinX*z33.XSinX+z31.SinX*z32.XCosX*z33.XCosX+z31.CosX*z32.XSinX*z33.XCosX+z31.CosX*z32.XCosX*z33.XSinX)...
    +((T11/T31)+(T13/T33))*(-z31.XSinX*z32.SinX*z33.XSinX+z31.XSinX*z32.CosX*z33.XCosX+z31.XCosX*z32.SinX*z33.XCosX+z31.XCosX*z32.CosX*z33.XSinX)

% calc E[(y*theta*cos(theta))]
E_ythC =  (T21/T31)*(z31.X2CosX*z32.CosX*z33.CosX-z31.X2CosX*z32.SinX*z33.SinX-z31.X2SinX*z32.CosX*z33.SinX-z31.X2SinX*z32.SinX*z33.CosX)...
    +(T22/T32)*(z31.CosX*z32.X2CosX*z33.CosX-z31.CosX*z32.X2SinX*z33.SinX-z31.SinX*z32.X2CosX*z33.SinX-z31.SinX*z32.X2SinX*z33.CosX)...
    +(T23/T33)*(z31.CosX*z32.CosX*z33.X2CosX-z31.CosX*z32.SinX*z33.X2SinX-z31.SinX*z32.CosX*z33.X2SinX-z31.SinX*z32.SinX*z33.X2CosX)...
    +((T21/T31)+(T22/T32))*(z31.XCosX*z32.XCosX*z33.CosX-z31.XCosX*z32.XSinX*z33.SinX-z31.XSinX*z32.XCosX*z33.SinX-z31.XSinX*z32.XSinX*z33.CosX)...
    +((T22/T32)+(T23/T33))*(z31.CosX*z32.XCosX*z33.XCosX-z31.CosX*z32.XSinX*z33.XSinX-z31.SinX*z32.XCosX*z33.XSinX-z31.SinX*z32.XSinX*z33.XCosX)...
    +((T21/T31)+(T23/T33))*(z31.XCosX*z32.CosX*z33.XCosX-z31.XCosX*z32.SinX*z33.XSinX-z31.XSinX*z32.CosX*z33.XSinX-z31.XSinX*z32.SinX*z33.XCosX)

% calc E[(y*theta*sin(theta))]
E_ythS =  (T21/T31)*(-z31.X2SinX*z32.SinX*z33.SinX+z31.X2SinX*z32.CosX*z33.CosX+z31.X2CosX*z32.SinX*z33.CosX+z31.X2CosX*z32.CosX*z33.SinX)...
    +(T22/T32)*(-z31.SinX*z32.X2SinX*z33.SinX+z31.SinX*z32.X2CosX*z33.CosX+z31.CosX*z32.X2SinX*z33.CosX+z31.CosX*z32.X2CosX*z33.SinX)...
    +(T23/T33)*(-z31.SinX*z32.SinX*z33.X2SinX+z31.SinX*z32.CosX*z33.X2CosX+z31.CosX*z32.SinX*z33.X2CosX+z31.CosX*z32.CosX*z33.X2SinX)...
    +((T21/T31)+(T22/T32))*(-z31.XSinX*z32.XSinX*z33.SinX+z31.XSinX*z32.XCosX*z33.CosX+z31.XCosX*z32.XSinX*z33.CosX+z31.XCosX*z32.XCosX*z33.SinX)...
    +((T22/T32)+(T23/T33))*(-z31.SinX*z32.XSinX*z33.XSinX+z31.SinX*z32.XCosX*z33.XCosX+z31.CosX*z32.XSinX*z33.XCosX+z31.CosX*z32.XCosX*z33.XSinX)...
    +((T21/T31)+(T23/T33))*(-z31.XSinX*z32.SinX*z33.XSinX+z31.XSinX*z32.CosX*z33.XCosX+z31.XCosX*z32.SinX*z33.XCosX+z31.XCosX*z32.CosX*z33.XSinX)```
コード（モンテカルロ法）
rng('default');
rng(1);
x = mvnrnd(mu_x,Sigma_x,1e8);

M_xC=mean(x(:,1).*cos(x(:,3)))
M_xS=mean(x(:,1).*sin(x(:,3)))
M_yC=mean(x(:,2).*cos(x(:,3)))
M_yS=mean(x(:,2).*sin(x(:,3)))
M_x2C=mean(x(:,1).*x(:,1).*cos(x(:,3)))
M_x2S=mean(x(:,1).*x(:,1).*sin(x(:,3)))
M_xyC=mean(x(:,1).*x(:,2).*cos(x(:,3)))
M_xyS=mean(x(:,1).*x(:,2).*sin(x(:,3)))
M_y2C=mean(x(:,2).*x(:,2).*cos(x(:,3)))
M_y2S=mean(x(:,2).*x(:,2).*sin(x(:,3)))
M_xthC=mean(x(:,1).*x(:,3).*cos(x(:,3)))
M_xthS=mean(x(:,1).*x(:,3).*sin(x(:,3)))
M_ythC=mean(x(:,2).*x(:,3).*cos(x(:,3)))
M_ythS=mean(x(:,2).*x(:,3).*sin(x(:,3)))
コード（結果の表示用）

e_xC=E_xC-M_xC
e_xS=E_xS-M_xS
e_yC=E_yC-M_yC
e_yS=E_yS-M_yS
e_xyC=E_xyC-M_xyC
e_xyS=E_xyS-M_xyS
e_x2C=E_x2C-M_x2C
e_x2S=E_x2S-M_x2S
e_y2C=E_y2C-M_y2C
e_y2S=E_y2S-M_y2S
e_xthC=E_xthC-M_xthC
e_xthS=E_xthS-M_xthS
e_ythC=E_ythC-M_ythC
e_ythS=E_ythS-M_ythS

sE = '\mathbb{E}'
s1 = sprintf(['| | %s[x\\\\cos(\\\\theta)] | %s[x\\\\sin(\\\\theta)] | %s[y\\\\cos(\\\\theta)] | %s[y\\\\sin(\\\\theta)] |' ...
' %s[xy\\\\cos(\\\\theta)] | %s[xy\\\\sin(\\\\theta)] | %s[x^2\\\\cos(\\\\theta)] | %s[x^2\\\\sin(\\\\theta)] | %s[y^2\\\\cos(\\\\theta)] | %s[y^2\\\\sin(\\\\theta)] |' ...
' %s[x\\\\theta\\\\cos(\\\\theta)] | %s[x\\\\theta\\\\sin(\\\\theta)] | %s[y\\\\theta\\\\cos(\\\\theta)] | %s[y\\\\theta\\\\sin(\\\\theta)] |\n'],...
sE,sE,sE,sE,sE,sE,sE,sE,sE,sE,sE,sE,sE,sE);
s2 = sprintf('| ---- | ---- | ---- | ---- | ---- | ---- | ---- | ---- | ---- | ---- | ---- | ---- | ---- | ---- | ---- |\n');
s3 = sprintf('|解析解|%.2f|%.2f|%.2f|%.2f|%.2f|%.2f|%.2f|%.2f|%.2f|%.2f|%.2f|%.2f|%.2f|%.2f|\n',...
E_xC,E_xS,E_yC,E_yS,E_xyC,E_xyS,E_x2C,E_x2S,E_y2C,E_y2S,E_xthC,E_xthS,E_ythC,E_ythS);
s4 = sprintf('|モンテカルロ|%.2f|%.2f|%.2f|%.2f|%.2f|%.2f|%.2f|%.2f|%.2f|%.2f|%.2f|%.2f|%.2f|%.2f|\n', ...
M_xC,M_xS,M_yC,M_yS,M_xyC,M_xyS,M_x2C,M_x2S,M_y2C,M_y2S,M_xthC,M_xthS,M_ythC,M_ythS);
s5 = sprintf('|差( \\\\times 10^{-4} )|%.2f|%.2f|%.2f|%.2f|%.2f|%.2f|%.2f|%.2f|%.2f|%.2f|%.2f|%.2f|%.2f|%.2f|\n', ...
e_xC1e4,e_xS1e4,e_yC1e4,e_yS1e4,e_xyC1e4,e_xyS1e4,e_x2C1e4,e_x2S1e4,e_y2C1e4,e_y2S1e4,e_xthC1e4,e_xthS1e4,e_ythC1e4,e_ythS1e4);
fprintf([s1,s2,s3,s4,s5])

おわりに

本記事で確率変数のモーメント計算はひと段落したと思う．
手計算は大変だが，一度計算しておくと正規分布の理解や変数の循環にも気づくためなかなか有用だったと思う．

これで準備はできたので，今後，Moment-based Kalman Filterについて取り組む．

Discussion