UA STAT675 統計計算I 隨機數生成2 線性遞歸模m

• • Multiple Recursive Generator (MRG)
  • MRG設計方法

線性遞歸模m是最常用的一種算法RNG(random number generator)，它的工作原理基礎是下面的遞歸式子
$$x_i=(a_1{x}_{i?1}+?+a_k{x}_{i?k})\mathrm{m}\mathrm{o}\mathrm{d}m$$

其中$m,k$均是自然數，$m$被稱為模數(modulus)，$k$被稱為階，$a_1,?,a_k\in {\mathbb{Z}}_m=\{0,1,?,m?1\}$，環${\mathbb{Z}}_m$上的運算在關于$m$同余的意義下封閉，當$m$是素數時，${\mathbb{Z}}_m=GL(m)$ (Galois域)。

上一講我們定義RNG的抽象結構：
$$RNG=(S,\mu ,f,U,g)$$

這里$S$是狀態集，$f$是狀態轉移函數，$f:S\to S$，$U$是輸出集，$U(0,1)$的RNG滿足$U=(0,1)$，$g:S\to U$是輸出函數，$f,g$都是確定性的函數，不產生隨機性，$\mu$是$S$的概率分布，是用來選擇初始狀態(initial state或者seed)的，記選擇的seed為$s_0$，則RNG的狀態轉移序列滿足：
$$s_t=f(s_{t?1})$$

從而輸出的隨機數為$\{u_t:u_t=g(u_t)\}$。

Multiple Recursive Generator (MRG)

Knuth (1998)證明了當$m$是素數時，最大的周期為$\rho =m^k?1$，當且僅當線性遞歸
$$x_i=(a_1{x}_{i?1}+?+a_k{x}_{i?k})\mathrm{m}\mathrm{o}\mathrm{d}m$$

的特征多項式
$$P(z)=z^k?a_1{z}^{k?1}??a_k$$

是Galois域$GL(m)$上的primitive polynomial；另一種等價敘述是
$$\underset{\nu \in {\mathbb{Z}}^{+}}{\inf }\{(z^{\nu }\mathrm{m}\mathrm{o}\mathrm{d}P(z))\mathrm{m}\mathrm{o}\mathrm{d}m=1\}=m^k?1$$

Knuth (1998)敘述了一種尋找primitive polynomial的簡化方法，如果$P(z)$是primitive polynomial，則遞歸的系數中除了$a_k$非零外，至少還有一個系數非零(記為$a_r$)；于是線性遞歸可以簡化為
$$x_n=(a_r{x}_{n?r}+a_k{x}_{n?k})\mathrm{m}\mathrm{o}\mathrm{d}m$$

即我們只需要確定$a_r,a_k$兩個系數的值，使得
$$P(z)=z^k?a_r{z}^{k?r}?a_k$$

是primitive polynomial即可。

根據這些結果，我們來構造隨機數生成器，為此我們需要確定狀態集、狀態轉移規則、輸出函數：

第$n$步的狀態集為
$$s_n=\{x_{n?k},x_{n?r},x_n\}$$

狀態轉移規則
$$x_n=(a_r{x}_{n?r}+a_k{x}_{n?k})\mathrm{m}\mathrm{o}\mathrm{d}m$$

輸出函數
$$u_n=\frac{x_n}{m}$$

這樣的隨機數生成器叫做Multiple Recursive Generator (MRG)。

一種特例是$k=1$時，

第$n$步的狀態集為
$$s_n=\{x_{n?1},x_n\}$$

狀態轉移規則
$$x_n=(a_1{x}_{n?1})\mathrm{m}\mathrm{o}\mathrm{d}m$$

輸出函數
$$u_n=\frac{x_n}{m}$$

這樣的隨機數生成器叫做linear congruential generator (LCG)。

MRG設計方法

要設計一個MRG，我們只需要確定系數$a_1,?,a_k,m$即可。注意到$x_n$與$m$都是整數，所以$u_n$是有理數，我們希望$u_n{～}_{iid}U(0,1)$，就需要$m$越大越好。另外，MRG涉及的主要計算就是在模$m$的意義下做線性遞歸，所以MRG的實現細節主要就是討論遞歸取模的計算。

Exact representation Method
第一種設計思路是讓遞歸式能夠被計算機準確表達，考慮在64-bit計算機IEEE float-point standard下的雙精度浮點數下進行模$m$計算，所有不超過$2^{53}$的整數都可以被準備表示，因此只要我們在設計參數時保證：
$$(|a_1|+?+|a_k|)(m?1)\le 2^{53}$$

$a_1{x}_{i?1}+?+a_k{x}_{i?k}$就可以被準確表示。

以遞歸中的一項為例，記$z_j=|a_j|x_{i?j}$，則$z_j\mathrm{m}\mathrm{o}\mathrm{d}m$用下面的算法計算：
$$\begin{gathered} y=|a_j|x_{i?j} \\ {z}_j=y?m?y/m? \end{gathered}$$

Approximate Factoring Method
以遞歸中的一項為例，記$z_j=|a_j|x_{i?j}$，取$|a_j|=i,i<\sqrt{m}$或者$|a_j|=?m/i?$，則$z_j\mathrm{m}\mathrm{o}\mathrm{d}m$用下面的算法計算：
$$\begin{gathered} q_j=?m/|a_j|? \\ {r}_j=m\mathrm{m}\mathrm{o}\mathrm{d}|a_j| \\ z=|a_j|(x_{i?j}?y{q}_j)?y{r}_j \end{gathered}$$

如果$z<0$，取$z_j=m+z$，否則$z_j=z$。

關于MRG的設計方法還有很多其他的觀點，感興趣的同學可以參考handbook of computational statistics chapter 3.

總結

以上是生活随笔為你收集整理的UA STAT675 统计计算I 随机数生成2 线性递归模m与Multiple Recursive Generator (MRG)的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。