5.5 可重現的亂數#

本節目標#

了解電腦亂數的真相：偽亂數是確定性演算法的輸出，seed 決定整條序列
認識 Mersenne Twister（Matsumoto & Nishimura, 1998）與 R 的 set.seed
說明為什麼科學計算需要可重現亂數（除錯、審查、回歸測試）
示範同一 seed 抽不同長度可以拼接：前 3 個 + 後 2 個 = 一次抽 5 個
用 rupy 把同一份 Python 編譯成 Rust 執行，驗證輸出與直譯完全相同

問題：你的模擬，別人跑得出同一份答案嗎？#

本章我們跑了大量模擬，每支程式第一行都是 g = r.set_seed(...)。如果拿掉這行，每次執行直方圖都會微妙地不同——作業交出去、論文投出去，審查者永遠無法重現你的圖。可重現性（reproducibility）不是潔癖，是科學計算的剛需：

除錯：模擬偶爾出錯、重跑又不錯？有 seed 才能把那次的亂數序列原樣叫回來。
論文審查：審查者用你的 seed 跑你的程式，應得到你的數字。
回歸測試：改了演算法之後跑同一 seed，輸出應該逐位元可比對，哪裡變了一目了然。

而「可重現」的根基是一個常被誤解的事實：

電腦沒有真亂數。set_seed(s) 之後的每一個「隨機」值都是確定性演算法算出來的——同一 seed 必得同一序列。這正是我們要的：偽亂數（pseudo-randomness）的缺點，恰好是科學計算需要的性質。

1998 — Mersenne Twister#

R 與 Python（random 模組）長年使用的預設亂數產生器是松本真與西村拓士 1998 年發表的 Mersenne Twister (MT19937)：週期長達 \(2^{19937}-1\) （恰好是梅森質數，故得名），輸出均勻性經過嚴格統計檢驗（diehard 測試）。它不可用於密碼學（看夠多輸出可反推內部狀態），但對模擬而言速度與品質俱佳，成為二十多年來統計軟體的事實標準。

R 的 set.seed(s) 用 seed 經固定程序初始化 MT19937 的內部狀態；rstatistics 的 r.set_seed(s) 對 R 的初始化程序逐位元等價——這意味著本書每個模擬都同時是 R 程式，換個殼能得到同一串數字。這在對照教科書（R 生態）與做可重現研究時價值巨大：同一段模擬，Python 與 R 的審查者都能驗證。

動手（一）：seed 的三個性質#

import rstatistics as r

# ---- (1) 同 seed 必得同序列（確定性）----
print("同 seed 兩次呼叫：")
print("  第一次:", r.set_seed(2026).rnorm(5))
print("  第二次:", r.set_seed(2026).rnorm(5))
print("  不同 seed:", r.set_seed(2027).rnorm(5))

# ---- (2) 拼接性：前 3 個 + 後 2 個 == 一次抽 5 個 ----
print()
a = r.set_seed(2026)
b = r.set_seed(2026)
print("一次抽 5 個        :", a.rnorm(5))
print("先抽 3 個再抽 2 個 :", b.rnorm(3), "+", b.rnorm(2))

# ---- (3) 位元等價：與 R 的 set.seed(42) 比對 ----
print()
print("r.set_seed(42) 的前兩個均勻亂數:", r.RNG(42).runif(2))
print("r.set_seed(42) 的前三個標準常態:", r.RNG(42).rnorm(3))

輸出（實際執行）：

同 seed 兩次呼叫：
  第一次: [0.5205890729185229, -1.07969076235228, 0.1392381150192734, -0.08474878494857653, -0.6666396152845961]
  第二次: [0.5205890729185229, -1.07969076235228, 0.1392381150192734, -0.08474878494857653, -0.6666396152845961]
  不同 seed: [-0.9348779593312257, -0.6938629367777, -0.5479450573071134, 0.7662709472389719, 1.8521001640647539]

一次抽 5 個        : [0.5205890729185229, -1.07969076235228, 0.1392381150192734, -0.08474878494857653, -0.6666396152845961]
先抽 3 個再抽 2 個 : [0.5205890729185229, -1.07969076235228, 0.1392381150192734] + [-0.08474878494857653, -0.6666396152845961]

r.set_seed(42) 的前兩個均勻亂數: [0.9148060434963554, 0.9370754132978618]
r.set_seed(42) 的前三個標準常態: [1.3709584471466685, -0.5646981713960888, 0.36312841133733914]

三件事：

確定性：同 seed 兩次，五個值完全相同。序列在 seed 設下的那一刻就全部決定了。
拼接性：rnorm(3) 之後接 rnorm(2)，等於一次 rnorm(5)。因為產生器是一台狀態機，每抽一個就推進一格——抽法怎麼分段不重要，消耗的亂數流量才重要。這也是為什麼模擬程式裡「誰先抽、抽幾個」的順序不能隨意改動：一改，後面所有數字全變。
位元等價：r.RNG(42).runif(2) 的值與 R 裡 set.seed(42); runif(2) 逐位元相同（可在 R 裡驗證，見下方對照）。本書每個模擬在 R 生态中都可以原樣重現。

動手（二）：rupy 編譯——同一份 .py，同一份輸出#

本書的 rstatistics 是純 Python；但如果模擬要跑幾億次，直譯會慢。工具 rupy（py2rs 專案）可以把同一份 Python 原始碼翻譯成 Rust 編譯執行。對可重現性而言，編譯前後輸出必須完全相同——我們真的來驗一次。

待編譯的程式 ttest.py：

import rstatistics as r

g = r.set_seed(42)
x = g.rnorm(30, mean=0.0, sd=1.0)
y = g.rnorm(25, mean=0.5, sd=1.0)
print(r.t_test(x, y))

直譯執行（CPython + 純 Python 的 rstatistics）：

cd /tmp && python3 ttest.py

輸出（實際執行）：

    Welch Two Sample t-test

data:  x and y
t = -1.1701, df = 52.75, p-value = 0.2472
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.9544103  0.2511616
sample estimates:
 mean of x  mean of y 
0.06858683 0.42021119

編譯成 Rust 執行（rupy 會把 rstatistics 連同 MT19937 亂數核心翻譯成 Rust 原始碼，再交給 cargo 編譯）：

cd /Users/Shared/ccc/project/rupy && ./target/release/rupy run /tmp/ttest.py

輸出（實際執行）：

    Welch Two Sample t-test

data:  x and y
t = -1.1701, df = 52.75, p-value = 0.2472
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.9544103  0.2511616
sample estimates:
 mean of x  mean of y 
0.06858683 0.42021119

逐字相同。想要看編譯出來的 Rust 長什麼樣，rupy emit 會印出（節錄）：

#[derive(Clone, Debug)]
pub struct rstatistics__RNG {
    pub mt: Vec<i64>,        // MT19937 的 624 個狀態字
    pub mti: i64,
    ...

MT19937 的內部狀態被原樣翻譯成 Rust struct——所以 seed 42 在 Rust 版推進出來的序列與 Python 版一致， \(t\) 檢定的每一步浮點運算也一致，輸出自然相同。rupy 另能編譯成 WebAssembly（rupy wasm），同一份原始碼跑在瀏覽器裡也是同一份輸出；本例用到較多執行期功能，此處僅以 Rust 為證。

這就是「可重現」的完整意涵：seed 定序列、序列定輸出、語言與平台只改速度不改答案。審查者拿到你的 .py 與 seed，用直譯、Rust、wasm 任何一條路，都會得到同一份 p 值。

與 R 對照#

驗證位元等價：在 R 裡跑同一個 seed，和 rstatistics 比對。

set.seed(42)
runif(2)
set.seed(42)
rnorm(3)

輸出（實際執行）：

[1] 0.9370754 0.2861395
[1]  1.3709584 -0.5646982  0.3631284

與 r.RNG(42).runif(2) = [0.9148060434963554, 0.9370754132978618]、r.RNG(42).rnorm(3) = [1.3709584471466685, -0.5646981713960888, 0.36312841133733914] 逐位元相同。

小結#

偽亂數是確定性演算法：seed 一設，整條序列就此注定——這是可重現性的來源。
亂數序列像一條帶子：抽法分段不重要，消耗流量才重要；改動抽樣順序會讓後面全部洗牌。
Mersenne Twister（Matsumoto & Nishimura, 1998）是 R / Python 的共同基礎；r.set_seed 與 R 位元等價，讓本書模擬可被 R 使用者直接驗證。
可重現性的最終考驗是跨實作：同一份 .py 經 rupy 編譯成 Rust（或 wasm），輸出與直譯逐字相同。
寫模擬程式的鐵律：第一行設 seed，並把它寫進論文/報告。

練習#

把 5.2 節 CLT 模擬的 set_seed(7) 拿掉，跑三次比較輸出；再放回去，確認每次執行結果一致。
用 r.RNG(1) 分別抽 rnorm(10)，以及先抽 runif(3) 再抽 rnorm(10)：兩者的 rnorm 序列相同嗎？為什麼？（提示：runif 與 rnorm 消耗同一條底層均勻序列，但消耗的量不同——rnorm 一個值常消耗 2 個均勻亂數。）
驗證拼接性對 rexp 也成立：g.rexp(5, 1.0) vs g.rexp(2, 1.0) + g.rexp(3, 1.0)。
把 5.3 節的蒙地卡羅積分程式（mc_integral）存成 .py，用 rupy run 編譯執行，比對直譯輸出。
找一個 seed，使 r.RNG(seed).rnorm(1)[0] > 2.5（用迴圈搜尋），並驗證該 seed 在 R 裡 set.seed(seed); rnorm(1) 得到同一個值。