hooked on mechatronics

いつまで経っても初心者

各確率分布に従う確率変数の4次までのモーメントとその計算方法およびプログラム

各確率分布に従う確率変数の4次までのモーメント

二項分布 B(n,p)

次数 モーメントE[Xn] 中心モーメントE[(X-E[X])n]
1次
2次
3次
4次

ポアソン分布

次数 モーメントE[Xn] 中心モーメントE[(X-E[X])n]
1次
2次
3次
4次

正規分布

次数 モーメントE[Xn] 中心モーメントE[(X-E[X])n]
1次
2次
3次
4次

ガンマ分布

次数 モーメントE[Xn] 中心モーメントE[(X-E[X])n]
1次
2次
3次
4次

n次モーメントの計算方法

 n 次モーメントは,確率変数  X積率母関数  M_X(t)=E[e^{tX}] n微分 t=0 を代入した値である(参考:積率母関数).
また, n 次の中心モーメントは

により, 1,\cdots,n 次モーメントから計算できる(二項定理).

プログラムによる計算

以上の計算方法を,MATLABのSymbolic Math Toolboxによって計算します.pythonのSymPyなどでも同様に計算できると思います.
以下はポアソン分布を例として用います.ポアソン分布の積率母関数  M_X(t)

であることを用いています.

他の分布についても,「変数の定義」と「積率母関数の定義」の部分を変更するだけで,モーメントを計算することができます.

%% 定義・初期化
% 変数の定義
t = sym('t');
la = sym('lambda');

% 求めるモーメントの次数の最大
n = 4;

% モーメントと中心モーメントを格納する変数の初期化
M_list = sym('M_list', [n, 1]);
central_M_list = sym('center_M_list', [n, 1]);

% 積率母関数の定義
syms M(t)
M(t) = exp(la*(exp(t)-1));

%% モーメント・中心モーメントの計算
for order = 1:n

    % order次モーメントの計算
    M_list(order) = subs(diff(M, order), t, 0);

    % order次の中心モーメントの計算
    if order == 1

        central_M_list(order) = M_list(1);

    else 

        % シグマ計算
        sum = 0;

        for k = 0 : order-1
    
            sum = sum + (-1)^k * nchoosek(order, k) * M_list(order - k) * M_list(1)^k;
    
        end

        sum = sum + (-1)^order * M_list(1)^order;
    
        % 結果を格納
        central_M_list(order) = sum;

    end

end

%% 結果の表示
disp("Moment")
disp(simplify(M_list))
disp("Central moment")
disp(simplify(central_M_list))