Saturday, July 29, 2006

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 twisterKISS 。但 Mersenne twister 較為普遍,可以很容易找到各種語言的版本,這裡強烈建議採用。

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

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

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

Testings of RNGs

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 分布

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

  • 機率一講機率一講

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

Tags: [] [] [] []

4 comments:

Yukuan said...

突然想到,JDK 1.1 以前附的 RNG 在一維的應用下良好,在二維的應用上會產生明顯的 pattern 。這個問題在 JDK 1.4 之後經實際測試已經被解掉了。兩次我都是為了要作一些模擬實驗才特別去著墨的。

Yukuan said...

順便一提, Python 的 random module 以 Mersenne Twister 法實作,當然是合格的。

Yukuan said...

後來查了 Python 下 Poisson Distribution 的現成 RNG 。發現 numarray.random_array 及 numpy.random 這兩個 module 都有提供。

其使用範例可以分別參考
Poisson Distribution (for a newbie)Numpy Example List

archiching said...

Very good introduction of randomness, especially it is writen in Chinese thanks very much.

I myself is studying architecture as well and researching the relationship between randomness and architecture. If you are interested we can keep contact and exchange idea.

Ching