8.2 隨機變數與 sympy.stats#

本節目標#

認識隨機變數、機率質量函數（pmf）、期望值與變異數
導覽 sympy.stats 模組：Die()、Coin()、E()、P()、variance()、density()
以「兩顆骰子之和」為例，對照手算與符號計算
實驗：模擬擲骰一萬次，畫直方圖並與理論值對照

隨機變數的概念#

隨機變數是把樣本空間的結果對應到數值的函數。例如擲兩顆骰子，隨機變數

\[ X = \text{兩顆點數之和} \]

把每個結果 \((i, j)\) 對應到 \(i + j\) 。離散隨機變數由其機率質量函數（pmf）完全描述：

\[ p(x) = P(X = x) \]

有了 pmf，兩個最重要的摘要量是期望值與變異數：

\[ E[X] = \sum_x x \, p(x), \qquad \mathrm{Var}(X) = E[(X - E[X])^2] = E[X^2] - (E[X])^2 \]

sympy.stats 模組導覽#

sympy.stats 把隨機變數變成 SymPy 的符號物件，之後所有計算（機率、期望值、變異數、密度）都是精確的符號運算。

from sympy.stats import Die, Coin, FiniteRV, E, P, variance, density
from sympy import Rational, Eq, simplify

# 骰子：預設 6 面
X = Die('X')
print(P(Eq(X, 3)))          # 1/6
print(P(X > 4))             # 1/3
print(E(X))                 # 7/2
print(variance(X))          # 35/12
print(density(X))           # {1: 1/6, 2: 1/6, ..., 6: 1/6}

# 偏斜的骰子：用 FiniteRV 指定各點數機率
from sympy import Rational
Y = FiniteRV('Y', {k: Rational(k, 21) for k in range(1, 7)})   # 加權骰子
print(density(Y))
print(E(Y).simplify())      # 13/3

# 銅板
C = Coin('C')
print(density(C))           # {H: 1/2, T: 1/2}
from sympy import Symbol
print(P(Eq(C, Symbol('H'))))    # 1/2

記得把 FiniteRV 加進 import：

from sympy.stats import FiniteRV

常用函數整理：

函數	用途
`P(條件)`	事件機率，如 `P(X > 3)`、`P(Eq(X, 2))`
`E(表達式)`	期望值，可作用在隨機變數的函數上，如 `E(X**2)`
`variance(X)`	變異數 \(\mathrm{Var}(X)\)
`std(X)`	標準差 \(\sigma\)
`cdf(X)(x)`	累積分布函數在 \(x\) 的值
`FiniteRV('Y', dict)`	自訂離散分布（自帶 pmf 字典）

範例：兩顆骰子之和的分布#

兩顆公平骰子之和 \(S = X_1 + X_2\) 的分布，手算是冗長的枚舉；SymPy 一行搞定：

from sympy.stats import Die, E, P, variance, density
from sympy import Rational

X1, X2 = Die('X1'), Die('X2')
S = X1 + X2                          # 隨機變數可以運算！

print(density(S))
# {2: 1/36, 3: 1/18, 4: 1/12, 5: 1/9, 6: 5/36, 7: 1/6,
#  8: 5/36, 9: 1/9, 10: 1/12, 11: 1/18, 12: 1/36}

print(P(Eq(S, 7)))       # 1/6，和為 7 的方法數最多（6 種）
print(P(S >= 10))        # 1/6
print(E(S))              # 7
print(variance(S))       # 35/6

注意 \(\mathrm{Var}(S) = 35/6 = 2 \times 35/12\) ，正好是單顆骰子變異數的兩倍——獨立隨機變數之和的變異數等於變異數之和（這是 8.4 節中央極限定理的伏筆）。

也可以驗證 \(E[S] = E[X_1] + E[X_2]\) 對任意（不必獨立的）隨機變數成立：

from sympy.stats import Die, E, variance
X1, X2 = Die('X1'), Die('X2')
print(E(X1 + X2) == E(X1) + E(X2))       # True
print(variance(X1 + X2) == variance(X1) + variance(X2))   # True（因為獨立）

期望值與變異數：手算 vs E() 對照#

例題：單顆骰子的期望值，手算為

\[ E[X] = \frac{1 + 2 + 3 + 4 + 5 + 6}{6} = \frac{7}{2} \]

變異數手算為

\[ E[X^2] = \frac{1 + 4 + 9 + 16 + 25 + 36}{6} = \frac{91}{6} \]

\[ \mathrm{Var}(X) = \frac{91}{6} - \left(\frac{7}{2}\right)^2 = \frac{91}{6} - \frac{49}{4} = \frac{35}{12} \]

用 SymPy 完整對照：

from sympy.stats import Die, E, variance, density
from sympy import Rational, nsimplify

X = Die('X')

# 手算版本
pmf = density(X)
E_manual  = sum(Rational(x) * p for x, p in pmf.items())
E2_manual = sum(Rational(x)**2 * p for x, p in pmf.items())
Var_manual = E2_manual - E_manual**2

print(E_manual)                 # 7/2
print(E2_manual)                # 91/6
print(Var_manual)               # 35/12

# SymPy 版本
print(E(X))                     # 7/2
print(E(X**2))                  # 91/6
print(variance(X))              # 35/12
print(variance(X) == Var_manual)   # True

# 定義式驗證：Var(X) = E[(X - E[X])^2]
from sympy.stats import sample
Y = Die('Y')
print(variance(Y))    # 35/12

E() 也能直接作用在隨機變數的任意函數上，例如 E((X - Rational(7, 2))**2) 同樣得 \(35/12\) ，這正是變異數的定義式。

動手實驗：擲骰一萬次的直方圖#

理論是精確的，實驗是頻率的近似。用 Python 內建的 random 模擬擲骰 10000 次，比較經驗頻率與理論機率 \(1/6\) ：

import random
from collections import Counter
from sympy import Rational

random.seed(42)                     # 固定亂數種子，結果可重現
N = 10_000

rolls = [random.randint(1, 6) for _ in range(N)]
counts = Counter(rolls)

print(f"{'點數':>4} {'次數':>6} {'頻率':>8} {'理論':>8}")
for face in range(1, 7):
    freq = counts[face] / N
    print(f"{face:>4} {counts[face]:>6} {freq:>8.4f} {Rational(1, 6).evalf(5):>8}")

再把「兩顆骰子之和」的模擬直方圖與 density() 的理論 pmf 並排比較：

import random
from collections import Counter
from sympy.stats import Die, density, E, variance
from sympy import Rational

random.seed(42)
N = 10_000
X1, X2 = Die('X1'), Die('X2')
S = X1 + X2

sums = [random.randint(1, 6) + random.randint(1, 6) for _ in range(N)]
counts = Counter(sums)

print(f"{'和':>3} {'模擬頻率':>9} {'理論 pmf':>9}")
for s in range(2, 13):
    freq = counts[s] / N
    theo = density(S)[s]
    print(f"{s:>3} {freq:>9.4f} {float(theo):>9.4f}")

print("模擬平均:", sum(sums) / N)        # ≈ 7.0
print("理論平均:", float(E(S)))          # 7.0
print("模擬變異數:", sum((x - 7)**2 for x in sums) / N)   # ≈ 5.83
print("理論變異數:", float(variance(S)))  # 5.8333...

若已安裝 matplotlib，可以畫出對比長條圖：

import random
import matplotlib.pyplot as plt
from collections import Counter
from sympy.stats import Die, density

random.seed(42)
S = Die('X1') + Die('X2')
sums = [random.randint(1, 6) + random.randint(1, 6) for _ in range(10_000)]
counts = Counter(sums)
theo = density(S)

s_range = range(2, 13)
plt.bar([s - 0.2 for s in s_range], [counts[s] / 10_000 for s in s_range],
        width=0.4, label='模擬')
plt.bar([s + 0.2 for s in s_range], [float(theo[s]) for s in s_range],
        width=0.4, label='理論')
plt.xticks(list(s_range))
plt.xlabel('兩顆骰子之和')
plt.ylabel('機率 / 頻率')
plt.legend()
plt.show()

試著修改實驗：

把 \(N\) 改成 \(100\) 、 \(1{,}000\) 、 \(100{,}000\) ，觀察模擬頻率如何隨樣本數收斂到理論值（大數法則）。
擲 5 顆骰子求和，理論期望值與變異數各是多少？用 E() 與 variance() 驗證，再模擬對照。
用自訂的偏斜骰子 FiniteRV('Z', {k: Rational(k, 21) for k in range(1, 7)}) 重做本節的對照實驗。

小結與練習#

隨機變數把結果對應到數值；離散隨機變數由 pmf 描述
sympy.stats 提供符號化的機率計算：P()、E()、variance()、density()、cdf()
兩顆骰子之和的分布可由 X1 + X2 直接得到，pmf 呈三角形、以 7 為峰
\(E[X_1 + X_2] = E[X_1] + E[X_2]\) 永遠成立；變異數相加則需要獨立性
模擬頻率隨樣本數增大而收斂到理論機率（大數法則的體驗）

練習

用 density() 求兩顆骰子「點數乘積」的期望值 \(E[X_1 X_2]\) ，並驗證 \(E[X_1 X_2] = E[X_1] E[X_2]\) （獨立性）。
擲三顆骰子，用 P() 求點數和為 10 的機率，與 8.1 節練習 1 的手算對照。
用 cdf() 計算 \(P(S \le 4)\) ，其中 \(S\) 是兩顆骰子之和，並手算驗證。
寫一個函數 roll_mean(n)，回傳擲 \(n\) 次骰子的平均；畫出 \(n = 1, 10, 100, 10000\) 時的平均值，觀察其向 \(3.5\) 收斂。
用 Python 的 random.randint(1, 6) 抽樣 100 次當作骰子樣本，計算樣本平均與樣本變異數，與理論值比較（除以 \(n\) 還是 \(n-1\) ？這是 8.4 節的主題）。