用Octave计算随机变量的矩
广告
{{v.name}}
矩是随机变量的一类基础数字特征，用来描述随机变量分布的形态特征（如集中趋势、离散程度、偏度、峰度等），分为原点矩和中心矩两大类，离散型和连续型随机变量的计算方式分别对应求和与积分。
原点矩的定义公式（以原点 \( x=0 \) 为中心）
设 \(X\) 为随机变量，\(k\) 为正整数，若 \(E(|X|^k) < \infty\)（保证期望存在），则称
\(\mu_k = E(X^k)\) 为 \(X\) 的 \(k\) 阶原点矩。
物理意义：类比质点系到原点的“矩”，\(k\) 阶原点矩是 \(X\) 取值的 \(k\) 次幂的加权平均。
低阶原点矩的熟悉形式：
- 1阶原点矩：\(\mu_1=E(X)\) → 就是数学期望，描述分布的集中趋势；
- 2阶原点矩：\(\mu_2=E(X^2)\) → 计算方差的关键中间量。
中心矩的定义公式（以数学期望 \( E(X) \) 为中心）
设 \(X\) 为随机变量，\(k\) 为正整数，若 \(E(|X|^k) < \infty\)（保证期望存在），则称
\(\mu_k' = E[(X - E(X))^k]\) 为 \(X\) 的 \(k\) 阶中心矩。
物理意义：类比质点系到质心的“矩”，\(k\) 阶中心矩是 \(X\) 偏离其数学期望的 \(k\) 次幂的加权平均。
低阶中心矩的熟悉形式：
- 1阶中心矩：\(\nu_1=E[X-E(X)]=E(X)-E(X)=0\) → 恒为0；
- 2阶中心矩：\(\nu_2=E[(X-E(X))^2]=D(X)\) → 就是方差，描述分布的离散程度；
- 3阶中心矩：\(\nu_3=E[(X-E(X))^3]\) → 用于计算偏度，描述分布的不对称性；
- 4阶中心矩：\(\nu_4=E[(X-E(X))^4]\) → 用于计算峰度，描述分布的陡峭程度。
离散型随机变量的矩计算
设离散型随机变量 \(X\) 的分布律为 \(P(X=x_i)=p_i\)（\(i=1,2,\dots\)），则：
\(k\) 阶原点矩：\(\mu_k = E(X^k) = \sum_{i=1}^{\infty} x_i^k p_i\)
\(k\) 阶中心矩：\(\nu_k = E\left[(X-E(X))^k\right] = \sum_{i=1}^{\infty} \left(x_i - E(X)\right)^k p_i\)
例如：0-1分布的矩计算
设 \(X\sim B(1,p)\)，分布律为 \(P(X=0)=1-p\)，\(P(X=1)=p\)。
- 1阶原点矩（期望）：\(\mu_1=0^1\cdot(1-p)+1^1\cdot p = p\)
- 2阶原点矩：\(\mu_2=0^2\cdot(1-p)+1^2\cdot p = p\)
- 2阶中心矩（方差）：\(\nu_2=(0-p)^2(1-p)+(1-p)^2 p = p(1-p)\)
- 3阶中心矩：\(\nu_3=(0-p)^3(1-p)+(1-p)^3 p = p(1-p)(1-2p)\)
程序代码如下
function ret = get_origin_moment_discrete(x_list, px_list, rank_value)
pkg load symbolic;
syms x;
n = length(x_list);
ret = 0;
for i = 1:n
ret = ret + (x_list(i) ^ rank_value) * px_list(i);
endfor
endfunction
function ret = get_central_moment_discrete(x_list, px_list, rank_value)
pkg load symbolic;
syms x;
n = length(x_list);
ret = 0;
EX = get_origin_moment_discrete(x_list, px_list, 1);
for i = 1:n
ret = ret + (x_list(i) - EX) ^ rank_value * px_list(i);
endfor
endfunction
>> p=sym('p');
>> get_origin_moment_discrete([0, 1], [1-p, p], 1)
ans = (sym) p
>> get_origin_moment_discrete([0, 1], [1-p, p], 2)
ans = (sym) p
>> get_central_moment_discrete([0, 1], [1-p, p], 2)
ans = (sym)
2 2
p ⋅(1 - p) + p⋅(1 - p)
>> get_central_moment_discrete([0, 1], [1-p, p], 3)
ans = (sym)
3 3
- p ⋅(1 - p) + p⋅(1 - p)连续型随机变量的矩计算
设连续型随机变量 \(X\) 的概率密度函数为 \(f(x)\)，则：
\(k\) 阶原点矩：\(\mu_k = E(X^k) = \int_{-\infty}^{+\infty} x^k f(x) dx\)
\(k\) 阶中心矩：\(\nu_k = E\left[(X-E(X))^k\right] = \int_{-\infty}^{+\infty} \left(x - E(X)\right)^k f(x) dx\)
例如：均匀分布 \(X\sim U(a,b)\) 的矩计算
\(X\) 的概率密度为 \(f(x)=\frac{1}{b-a}\)（\(a\le x\le b\)）。
- 1阶原点矩（期望）：\(\mu_1=\int_a^b x\cdot\frac{1}{b-a}dx = \frac{a+b}{2}\)
- 2阶原点矩：\(\mu_2=\int_a^b x^2\cdot\frac{1}{b-a}dx = \frac{a^2+ab+b^2}{3}\)
- 2阶中心矩（方差）：\(\nu_2=\int_a^b \left(x-\frac{a+b}{2}\right)^2\cdot\frac{1}{b-a}dx = \frac{(b-a)^2}{12}\)
程序代码如下
function ret = get_origin_moment_continuous(f, x, rank_value, a, b)
pkg load symbolic;
try
ret = double(int((x ^ rank_value) * f, x, a, b));
catch
ret = int((x ^ rank_value) * f, x, a, b);
end_try_catch
endfunction
function ret = get_central_moment_continuous(f, x, rank_value, a, b)
pkg load symbolic;
EX = get_origin_moment_continuous(f, x, 1, a, b);
try
ret = double(int(((x - EX) ^ rank_value) * f, x, a, b));
catch
ret = int(((x - EX) ^ rank_value) * f, x, a, b);
end_try_catch
endfunction
>> a = sym('a');
>> b = sym('b');
>> x = sym('x');
>> get_origin_moment_continuous((1 / (b-a)), x, 1, a, b)
ans = (sym)
2 2
a b
-------- - --------
2⋅a - 2⋅b 2⋅a - 2⋅b
>> get_origin_moment_continuous((1 / (b-a)), x, 2, a, b)
ans = (sym)
3 3
a b
-------- - --------
3⋅a - 3⋅b 3⋅a - 3⋅b
>> get_central_moment_continuous((1 / (b-a)), x, 2, a, b)
ans = (sym)
3 2 ⎛ 2 2⎞ 3 2 ⎛ 2 2⎞
a a ⋅(a + b) a⋅⎝- a - 2⋅a⋅b - b ⎠ b b ⋅(a + b) b⋅⎝- a - 2⋅a⋅b - b ⎠
-------- - --------- - --------------------- - -------- + -------- + ---------------------------
3⋅a - 3⋅b 2⋅a - 2⋅b 4⋅a - 4⋅b 3⋅a - 3⋅b 2⋅a - 2⋅b 4⋅a - 4⋅b
