《非一様乱数》のソースを表示

'''【ひいちようらんすう (non-uniform random numbers) 】'''

　一様分布以外の確率分布に従う乱数を[[非一様乱数]]という. これらの乱数は, 一様乱数に適当な変換を施して作るのがふつうである. 変換方法には, 種々の分布に対して適用可能な一般的なものと, 個々の分布の特徴を利用するものとがある. これらの方法を網羅的に扱っているのが [3]である. ここでは, 一般的な変換法を二つ紹介した後, 実用上大事な非一様乱数生成法を述べる. 以下では, (整数型の)一様乱数を正規化して得られる区間[0,1)上の(実数型の)一様乱数を<math>U, U_1, U_2\, </math>等で表し, これらは互いに独立であるものと仮定する. また, 発生させたい乱数を<math>X\, </math>とし, その分布関数を<math>F(x)\, </math>, それが連続分布ならば, その密度関数を<math>f(x)\, </math>で表す. 

[逆関数法]

　<math>F(x)\, </math>の逆関数<math>F^{-1}(x)\, </math>を使って<math>X=F^{-1}(U)\, </math>とするものである. 例えば, 指数分布, ワイブル分布, ロジスティック分布(スケール・パラメータはいずれも1とする)なら, それぞれ<math>X=-\log U\, </math>, <math>X=(-\log U)^\alpha\, </math>, <math>X=\log(U/(1-U))\, </math>とすればよい. 逆関数が解析的に求まらない場合には, 近似式等を用いることになる. 

[二者択一法]

　ウォーカー(A. J. Walker) [4] によって提案されたもので, 別名法(alias method)ともいう. 

　取りうる値の個数<math>(n)\, </math>が有限個の離散分布に従う乱数を, <math>n\, </math>に無関係なO(1)の速度で生成できる, 簡単で効率的な方法である. ただし, 初期設定にO<math>(n)\, </math>の時間とO<math>(n)\, </math>のメモリを必要とする. 例えばポアソン分布のように, 取りうる値の個数が無限の場合には, 分布の裾を適当に打ち切って適用すればよい. 

[正規分布]

　逆関数法を適用する場合には, 逆関数<math>F^{-1}(x)\, </math>が解析的には求まらないので, 例えば次の山内二郎の近似式を使う. 


<center>
<math>
\begin{array}{rll}
F^{-1}(x)& = &\left\{
 \begin{array}{ll}
  -w, & \ \ \ \ \  0 < x < 0.5 \\ 
  +w, & \ \ \ \ \  0.5 \leq x < 1.0
 \end{array}\right.\\ 
 w & = & \{ z(2.0611786-\frac{5.7262204}{z+11.640595})\}^{1/2}\\ 
 z & = & -\log\{4x(1-x)\}
\end{array}
\, </math>
</center>


ボックス・ミュラー(Box-Muller)法は, 変換公式


<center>
<math>
X_1=\sqrt{-2\log U_1}\cos(2\pi U_2),
\, </math>
</center>

<center>
<math>
X_2=\sqrt{-2\log U_1}\sin(2\pi U_2)
\, </math>
</center>


を使って, 2個の一様乱数<math>U_1,U_2\, </math>から互いに独立な標準正規乱数<math>X_1,X_2\, </math>を作るものである. 

　平均ベクトルが<math>\boldsymbol{\mu}\, </math>, 分散共分散行列が<math>V\, </math>の<math>n\, </math>次元正規乱数ベクトル<math>{\mathbf{ Y}}\, </math>を生成するためには, 互いに独立な<math>n\, </math>個の標準正規乱数からなるベクトル<math>{\mathbf{ X}}\, </math>に変換


<center>
<math>
\mathbf{Y} = \boldsymbol{\mu} + A \mathbf{X}
\, </math>
</center>


を施せばよい. ただし, ベクトルはすべて列ベクトルとし, <math>A\, </math>は<math>V=AA^\top\, </math>(<math>A^\top\, </math>は<math>A\, </math>の転置行列)を満たす下三角行列であり, コレスキー(Cholesky)分解によって求める. 

[アーラン分布]

　フェイズが<math>k\, </math>(正整数)のアーラン分布の密度関数は


<center>
:<math>
f(x)=\lambda \mathrm{e}^{-\lambda x}(\lambda x)^{k-1}/(k-1)!
\ \ \ \ \ \ \ (x\geq 0)
\, </math>
</center>


である. これは, パラメータ(平均値の逆数)が<math>\lambda\, </math>の指数分布に従う<math>k\, </math>個の確率変数の和の分布であるから, 


<center>
<math>
X=-\log(U_1\cdots U_k)/\lambda
\, </math>
</center>


とすればよい. 

[ポアソン分布]

　平均値が<math>\lambda\, </math>のポアソン乱数は, 平均が<math>1/\lambda\, </math>の指数乱数との関係を利用して作ることができる. すなわち, <math>U_1U_2\cdots U_n > {\mbox{e}}^{-\lambda}\, </math>が成り立つ最大の整数<math>n\, </math>を<math>X\, </math>とすればよい. 1個のポアソン乱数を発生するのに必要な一様乱数の個数の平均値は<math>\lambda+1\, </math>であるから, <math>\lambda\, </math>が大変に大きいときには, この方法は効率が悪いので, 二者択一法を使うほうがよい. 



----
'''参考文献'''

[1] 伏見正則, 『乱数』(UP応用数学選書12), 東京大学出版会, 1989. 

[2] D. E. Knuth, ''The Art of Computer Programming, Vol.2: Seminumerical Algorithms, 2nd ed.,'' Addison-Wesley, 1981. 渋谷政昭訳, 『準数値算法/乱数』, サイエンス社, 1981.

[3] L. Devroye, ''Non-Uniform Random Variate Generation'', Springer-Verlag, 1986.

[4] A. J. Walker, "An Efficient Method for Generating Discrete Random Variables with General Distributions," ''ACM Transactions on Mathematical Software'', '''3''' (1977), 253-256.