Not to be Fooled by Randomness

最近抱讀了一些機率的書，突然發現這門學問還滿有趣的。除了有趣，一些基本的機率分配及常見的隨機過程都非常實用。在複習的過程，當然也要作個紀錄：

Random Number Generators

如果想跑些像 Monte Carlo 這類對隨機性很要求的 simulation ，有必要好好檢視一下手邊的隨機數產生器（Random Number Generator, RNG）是否合格。

線性同餘法（Linear congruential generators, LCGs）是多數開發環境標準的 RNG。它雖可以輕快、方便地產生隨機數；不幸地，它產生的隨機數，很容易就不夠隨機，不適合用在嚴肅的場合。

對於 RNG ，我們重視的有：

1. 演算法要簡單、計算要快速
2. 要能通過各種隨機性測試，特定模式不可出現在我們運用的場合。
3. 用於密碼學領域時，還要夠安全，不易被破解

最早被採用的平方取中（middle-squre method）無法通過大部份的隨機性檢驗。廣為採用的線性同餘法（Linear congruential generators）很容易就因為參數沒選對，而產生差勁的隨機序列。

在隨機性要求高的場合，可以採用 Mersenne twister 或 KISS 。但 Mersenne twister 較為普遍，可以很容易找到各種語言的版本，這裡強烈建議採用。

有一種叫 Linear feedback shift register, LFSR 的，雖然可以方便用軟體撰寫，但其設計是為了可以用簡單的數位電路實作，以快速地產生隨機數。

另一方面，在密碼學這類還要求安全性的場合，Blum-Blum-Shub, BBS 是一類很經典的方法，此外還有 Fortuna, Yarrow, 及 ISAAC 等方法可用。

還有一種叫做 HAVEGE 的，安全性也很高，它會即時讀取電腦內如 cache, pipeline states 或 TLB 這類快數轉變的資料，缺點是不好驗證及除錯。

Testings of RNGs

• The NIST Collection
• Diehard Tests by George Marsaglia

Conversions

上述的 RNG 產生的隨機數都是落在 [0, n] 的均勻（uniform）整數（n 表示隨機數的上限）；它用來虛擬隨機變數（random variable） X_{_int}，亦即：

X_{_int} -> Unif[0, n]

其 p.m.f 如下：

p(x)= 1/n, if 0 <= x <= n, x in N;
p(x)= 0, otherwise

想把它轉換成 [0, 1) 的 uniform 實數，可以很簡單地以下式達成：

X_{_real}: [0, n] in N -> [0, 1) in R
X_{_real}(X_{_int})= X_{_int} / (n+1.0)
X_{_real} -> Unif[0, 1)

其 p.d.f：

f(x)= 1, if 0 <= x < 1;
f(x)= 0, otherwise

要轉到 [m, n] 的 uniform 整數也很簡單：

X_{_IntRng}: [0, 1) in R -> [m, n] in N
X_{_IntRng}(X_{_real})= m + floor((n-m+1) * X_{_real})
X_{_IntRng} -> Unif[m, n]

其 p.m.f：

p(x)= 1/(n-m), if m <= x <= n, x in N;
p(x)= 0, otherwise

想將分佈由 uniform 轉換成 non-uniform ，基本概念就是：

1. 以 p.d.f 或 p.m.f 求其 c.d.f ，也就是由 f(x) 或 p(x) 算出 F(x)
2. 求出 c.d.f 的反函數 F^-1(y)
3. 利用 uniform 的 RNG 取得 y
4. 把 y 餵給 F^-1(y) 就得到一個 non-uniform random number

如果我們要的 non-uniform random number 是離散的（discrete），還可以用 alias method 或 table lookup methods 這些更快的方法。

如果我們要的是 continuous non-uniform random number，很可能遇到反函數 F^-1(y) 不好求得的問題。這時可以 acceptance-rejection method 來抽樣，缺點是如果 rejection 發生的頻率太高，會嚴重影響執行效率。所以實用上大都採用拼湊法（Monty Python）及其變形。

在實用上，有沒有較省事的方案呢？如果開發環境是 C++ 這裡強烈建議使用 Boost Random Number Library，它很快就會被納入 C++ Standard Library 了。如果使用的是 C ，那 GNU Scientific Library 附的 Random Number Generation 是不錯的選擇。

Book List

講機率的書非常多，這裡只列出兩本小品文：

• 大於1/2──投資、愛情、生活的獲勝機率

  

  怎樣可以稱為隨機？公車為什麼總是來得比平均時間慢？生日相同的機率有多大？如何戀愛成功？賭徒的命運為何總是步向毀滅？隨機漫步看起來竟然一點都不隨機？……一連串的問題，本書都作了探討。

• 隨機的致富陷阱

  

  這本書是幾年前由 Claude 推薦的。作者強調真正隨機的現象看起來很可能一點都不隨機；而看起來明顯不隨機的現象，很可能是隨機下的必然結果。許多領域，如醫藥界，其隨機分佈充滿偏態及不對稱，所以平均數或中位數無法傳達多少有意義的事。稀有事件總是會發生，但統計學家總是察覺不到。更慘的是，許多現象，其隨機分佈不是定常的（stationary），它會隨外界的事件演化，這時標準的統計方法更是失靈了。

Cool Links

• 二項分布與大數法則

  上面這種機率分布稱為二項分布。一般的二項分布是這樣的： 假設某事件的發生率為 p，而試驗做了 n 次。這 n 次中，某事件發生 x 次的機率為

• 機率與訊息

  本文的目的有二：一來淺介訊息理論，二來藉訊息理論之介紹說明機率論的方法。 首先，我們回憶一下機率的基本概念。

• 漫談布朗運動

  西元1827年，英國植物學家勞伯‧布朗 (Robert Brown) 利用一般的顯微鏡觀察懸浮於水中的花粉粒時，發現這些花粉粒會做連續快速而不規則的隨機移動，這種移動稱為布朗運動 (Brownian motion)。

• 馬可夫鏈的簡介

  「不是互相獨立」也就是說互相關聯的意思，但是要樣相關呢？如何在相關中作一些簡單的分類呢？馬可夫鏈就是要描述在「相關」這個概念中最簡單的一種。

• Poisson 分布

  二項分布是離散型機率模型中最有名的一個，其次是 Poisson 分布，它可以看成為二項分布的一種極限情形。

• 機率一講

  機率論的一個發源地是記述統計學；明瞭記述統計學的概念，對於機率論的學習，大有助益！

Tags: [math] [random] [probability] [RNG]