「《非一様乱数》」の版間の差分
1行目: | 1行目: | ||
'''【ひいちようらんすう (non-uniform random numbers) 】''' | '''【ひいちようらんすう (non-uniform random numbers) 】''' | ||
− | 一様分布以外の確率分布に従う乱数を[[非一様乱数]]という. これらの乱数は, 一様乱数に適当な変換を施して作るのがふつうである. 変換方法には, 種々の分布に対して適用可能な一般的なものと, 個々の分布の特徴を利用するものとがある. これらの方法を網羅的に扱っているのが [3]である. ここでは, 一般的な変換法を二つ紹介した後, 実用上大事な非一様乱数生成法を述べる. 以下では, (整数型の)一様乱数を正規化して得られる区間[0,1)上の(実数型の)一様乱数を$U, U_1, U_2$等で表し, これらは互いに独立であるものと仮定する. また, 発生させたい乱数を$X$とし, その分布関数を$F(x)$, それが連続分布ならば, その密度関数を$f(x)$で表す. | + | 一様分布以外の確率分布に従う乱数を[[非一様乱数]]という. これらの乱数は, 一様乱数に適当な変換を施して作るのがふつうである. 変換方法には, 種々の分布に対して適用可能な一般的なものと, 個々の分布の特徴を利用するものとがある. これらの方法を網羅的に扱っているのが [3]である. ここでは, 一般的な変換法を二つ紹介した後, 実用上大事な非一様乱数生成法を述べる. 以下では, (整数型の)一様乱数を正規化して得られる区間[0,1)上の(実数型の)一様乱数を$<math>U, U_1, U_2</math>$等で表し, これらは互いに独立であるものと仮定する. また, 発生させたい乱数を$<math>X</math>$とし, その分布関数を$<math>F(x)</math>$, それが連続分布ならば, その密度関数を$<math>f(x)</math>$で表す. |
[逆関数法] | [逆関数法] | ||
− | $F(x)$の逆関数$F^{-1}(x)$を使って$X=F^{-1}(U)$とするものである. 例えば, 指数分布, ワイブル分布, ロジスティック分布(スケール・パラメータはいずれも1とする)なら, それぞれ$X=-\log U$, $X=(-\log U)^\alpha$, $X=\log(U/(1-U))$とすればよい. 逆関数が解析的に求まらない場合には, 近似式等を用いることになる. | + | $<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>$とすればよい. 逆関数が解析的に求まらない場合には, 近似式等を用いることになる. |
[二者択一法] | [二者択一法] | ||
11行目: | 11行目: | ||
ウォーカー(A. J. Walker) [4] によって提案されたもので, 別名法(alias method)ともいう. | ウォーカー(A. J. Walker) [4] によって提案されたもので, 別名法(alias method)ともいう. | ||
− | 取りうる値の個数$(n)$が有限個の離散分布に従う乱数を, $n$に無関係なO(1)の速度で生成できる, 簡単で効率的な方法である. ただし, 初期設定にO$(n)$の時間とO$(n)$のメモリを必要とする. 例えばポアソン分布のように, 取りうる値の個数が無限の場合には, 分布の裾を適当に打ち切って適用すればよい. | + | 取りうる値の個数$<math>(n)</math>$が有限個の離散分布に従う乱数を, $<math>n</math>$に無関係なO(1)の速度で生成できる, 簡単で効率的な方法である. ただし, 初期設定にO$<math>(n)</math>$の時間とO$<math>(n)</math>$のメモリを必要とする. 例えばポアソン分布のように, 取りうる値の個数が無限の場合には, 分布の裾を適当に打ち切って適用すればよい. |
[正規分布] | [正規分布] | ||
− | 逆関数法を適用する場合には, 逆関数$F^{-1}(x)$が解析的には求まらないので, 例えば次の山内二郎の近似式を使う. | + | 逆関数法を適用する場合には, 逆関数$<math>F^{-1}(x)</math>$が解析的には求まらないので, 例えば次の山内二郎の近似式を使う. |
45行目: | 45行目: | ||
を使って, 2個の一様乱数$U_1,U_2$から互いに独立な標準正規乱数$X_1,X_2$を作るものである. | を使って, 2個の一様乱数$U_1,U_2$から互いに独立な標準正規乱数$X_1,X_2$を作るものである. | ||
− | 平均ベクトルが$\mbox{\boldsymbol$\mu$}$, 分散共分散行列が$V$の$n$次元正規乱数ベクトル${\mathbf{ Y}}$を生成するためには, 互いに独立な$n$個の標準正規乱数からなるベクトル${\mathbf{ X}}$に変換 | + | 平均ベクトルが$<math>\mbox{\boldsymbol$\mu$}</math>$, 分散共分散行列が$<math>V</math>$の$<math>n</math>$次元正規乱数ベクトル$<math>{\mathbf{ Y}}</math>$を生成するためには, 互いに独立な$<math>n</math>$個の標準正規乱数からなるベクトル$<math>{\mathbf{ X}}</math>$に変換 |
53行目: | 53行目: | ||
− | を施せばよい. ただし, ベクトルはすべて列ベクトルとし, $A$は$V=AA^\top$($A^\top$は$A$の転置行列)を満たす下三角行列であり, コレスキー(Cholesky)分解によって求める. | + | を施せばよい. ただし, ベクトルはすべて列ベクトルとし, $<math>A</math>$は$<math>V=AA^\top</math>$($<math>A^\top</math>$は$<math>A</math>$の転置行列)を満たす下三角行列であり, コレスキー(Cholesky)分解によって求める. |
[アーラン分布] | [アーラン分布] | ||
− | フェイズが$k$(正整数)のアーラン分布の密度関数は | + | フェイズが$<math>k</math>$(正整数)のアーラン分布の密度関数は |
66行目: | 66行目: | ||
− | である. これは, パラメータ(平均値の逆数)が$\lambda$の指数分布に従う$k$個の確率変数の和の分布であるから, | + | である. これは, パラメータ(平均値の逆数)が$<math>\lambda</math>$の指数分布に従う$<math>k</math>$個の確率変数の和の分布であるから, |
78行目: | 78行目: | ||
[ポアソン分布] | [ポアソン分布] | ||
− | 平均値が$\lambda$のポアソン乱数は, 平均が$1/\lambda$の指数乱数との関係を利用して作ることができる. すなわち, $U_1U_2\cdots U_n>{\mbox{\rm e}}^{-\lambda}$が成り立つ最大の整数$n$を$X$とすればよい. 1個のポアソン乱数を発生するのに必要な一様乱数の個数の平均値は$\lambda+1$であるから, $\lambda$が大変に大きいときには, この方法は効率が悪いので, 二者択一法を使うほうがよい. | + | 平均値が$<math>\lambda</math>$のポアソン乱数は, 平均が$<math>1/\lambda</math>$の指数乱数との関係を利用して作ることができる. すなわち, $<math>U_1U_2\cdots U_n>{\mbox{\rm e}}^{-\lambda}</math>$が成り立つ最大の整数$<math>n</math>$を$<math>X</math>$とすればよい. 1個のポアソン乱数を発生するのに必要な一様乱数の個数の平均値は$<math>\lambda+1</math>$であるから, $<math>\lambda</math>$が大変に大きいときには, この方法は効率が悪いので, 二者択一法を使うほうがよい. |
2007年7月12日 (木) 13:22時点における版
【ひいちようらんすう (non-uniform random numbers) 】
一様分布以外の確率分布に従う乱数を非一様乱数という. これらの乱数は, 一様乱数に適当な変換を施して作るのがふつうである. 変換方法には, 種々の分布に対して適用可能な一般的なものと, 個々の分布の特徴を利用するものとがある. これらの方法を網羅的に扱っているのが [3]である. ここでは, 一般的な変換法を二つ紹介した後, 実用上大事な非一様乱数生成法を述べる. 以下では, (整数型の)一様乱数を正規化して得られる区間[0,1)上の(実数型の)一様乱数を$構文解析に失敗 (MathML、ただし動作しない場合はSVGかPNGで代替(最新ブラウザーや補助ツールに推奨): サーバー「https://en.wikipedia.org/api/rest_v1/」から無効な応答 ("Math extension cannot connect to Restbase."):): {\displaystyle U, U_1, U_2} $等で表し, これらは互いに独立であるものと仮定する. また, 発生させたい乱数を$構文解析に失敗 (MathML、ただし動作しない場合はSVGかPNGで代替(最新ブラウザーや補助ツールに推奨): サーバー「https://en.wikipedia.org/api/rest_v1/」から無効な応答 ("Math extension cannot connect to Restbase."):): {\displaystyle X} $とし, その分布関数を$$, それが連続分布ならば, その密度関数を$構文解析に失敗 (MathML、ただし動作しない場合はSVGかPNGで代替(最新ブラウザーや補助ツールに推奨): サーバー「https://en.wikipedia.org/api/rest_v1/」から無効な応答 ("Math extension cannot connect to Restbase."):): {\displaystyle f(x)} $で表す.
[逆関数法]
$$の逆関数$構文解析に失敗 (MathML、ただし動作しない場合はSVGかPNGで代替(最新ブラウザーや補助ツールに推奨): サーバー「https://en.wikipedia.org/api/rest_v1/」から無効な応答 ("Math extension cannot connect to Restbase."):): {\displaystyle F^{-1}(x)} $を使って$$とするものである. 例えば, 指数分布, ワイブル分布, ロジスティック分布(スケール・パラメータはいずれも1とする)なら, それぞれ$構文解析に失敗 (MathML、ただし動作しない場合はSVGかPNGで代替(最新ブラウザーや補助ツールに推奨): サーバー「https://en.wikipedia.org/api/rest_v1/」から無効な応答 ("Math extension cannot connect to Restbase."):): {\displaystyle X=-\log U} $, $$, $構文解析に失敗 (MathML、ただし動作しない場合はSVGかPNGで代替(最新ブラウザーや補助ツールに推奨): サーバー「https://en.wikipedia.org/api/rest_v1/」から無効な応答 ("Math extension cannot connect to Restbase."):): {\displaystyle X=\log(U/(1-U))} $とすればよい. 逆関数が解析的に求まらない場合には, 近似式等を用いることになる.
[二者択一法]
ウォーカー(A. J. Walker) [4] によって提案されたもので, 別名法(alias method)ともいう.
取りうる値の個数$構文解析に失敗 (MathML、ただし動作しない場合はSVGかPNGで代替(最新ブラウザーや補助ツールに推奨): サーバー「https://en.wikipedia.org/api/rest_v1/」から無効な応答 ("Math extension cannot connect to Restbase."):): {\displaystyle (n)} $が有限個の離散分布に従う乱数を, $構文解析に失敗 (MathML、ただし動作しない場合はSVGかPNGで代替(最新ブラウザーや補助ツールに推奨): サーバー「https://en.wikipedia.org/api/rest_v1/」から無効な応答 ("Math extension cannot connect to Restbase."):): {\displaystyle n} $に無関係なO(1)の速度で生成できる, 簡単で効率的な方法である. ただし, 初期設定にO$構文解析に失敗 (MathML、ただし動作しない場合はSVGかPNGで代替(最新ブラウザーや補助ツールに推奨): サーバー「https://en.wikipedia.org/api/rest_v1/」から無効な応答 ("Math extension cannot connect to Restbase."):): {\displaystyle (n)} $の時間とO$構文解析に失敗 (MathML、ただし動作しない場合はSVGかPNGで代替(最新ブラウザーや補助ツールに推奨): サーバー「https://en.wikipedia.org/api/rest_v1/」から無効な応答 ("Math extension cannot connect to Restbase."):): {\displaystyle (n)} $のメモリを必要とする. 例えばポアソン分布のように, 取りうる値の個数が無限の場合には, 分布の裾を適当に打ち切って適用すればよい.
[正規分布]
逆関数法を適用する場合には, 逆関数$構文解析に失敗 (MathML、ただし動作しない場合はSVGかPNGで代替(最新ブラウザーや補助ツールに推奨): サーバー「https://en.wikipedia.org/api/rest_v1/」から無効な応答 ("Math extension cannot connect to Restbase."):): {\displaystyle F^{-1}(x)} $が解析的には求まらないので, 例えば次の山内二郎の近似式を使う.
- 構文解析に失敗 (MathML、ただし動作しない場合はSVGかPNGで代替(最新ブラウザーや補助ツールに推奨): サーバー「https://en.wikipedia.org/api/rest_v1/」から無効な応答 ("Math extension cannot connect to Restbase."):): {\displaystyle \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} }
ボックス・ミュラー(Box-Muller)法は, 変換公式
- 構文解析に失敗 (MathML、ただし動作しない場合はSVGかPNGで代替(最新ブラウザーや補助ツールに推奨): サーバー「https://en.wikipedia.org/api/rest_v1/」から無効な応答 ("Math extension cannot connect to Restbase."):): {\displaystyle X_1=\sqrt{-2\log U_1}\cos(2\pi U_2), }
- 構文解析に失敗 (MathML、ただし動作しない場合はSVGかPNGで代替(最新ブラウザーや補助ツールに推奨): サーバー「https://en.wikipedia.org/api/rest_v1/」から無効な応答 ("Math extension cannot connect to Restbase."):): {\displaystyle X_2=\sqrt{-2\log U_1}\sin(2\pi U_2) }
を使って, 2個の一様乱数$U_1,U_2$から互いに独立な標準正規乱数$X_1,X_2$を作るものである.
平均ベクトルが$構文解析に失敗 (構文エラー): {\displaystyle \mbox{\boldsymbol$\mu$}} $, 分散共分散行列が$構文解析に失敗 (MathML、ただし動作しない場合はSVGかPNGで代替(最新ブラウザーや補助ツールに推奨): サーバー「https://en.wikipedia.org/api/rest_v1/」から無効な応答 ("Math extension cannot connect to Restbase."):): {\displaystyle V} $の$構文解析に失敗 (MathML、ただし動作しない場合はSVGかPNGで代替(最新ブラウザーや補助ツールに推奨): サーバー「https://en.wikipedia.org/api/rest_v1/」から無効な応答 ("Math extension cannot connect to Restbase."):): {\displaystyle n} $次元正規乱数ベクトル$構文解析に失敗 (MathML、ただし動作しない場合はSVGかPNGで代替(最新ブラウザーや補助ツールに推奨): サーバー「https://en.wikipedia.org/api/rest_v1/」から無効な応答 ("Math extension cannot connect to Restbase."):): {\displaystyle {\mathbf{ Y}}} $を生成するためには, 互いに独立な$構文解析に失敗 (MathML、ただし動作しない場合はSVGかPNGで代替(最新ブラウザーや補助ツールに推奨): サーバー「https://en.wikipedia.org/api/rest_v1/」から無効な応答 ("Math extension cannot connect to Restbase."):): {\displaystyle n} $個の標準正規乱数からなるベクトル$構文解析に失敗 (MathML、ただし動作しない場合はSVGかPNGで代替(最新ブラウザーや補助ツールに推奨): サーバー「https://en.wikipedia.org/api/rest_v1/」から無効な応答 ("Math extension cannot connect to Restbase."):): {\displaystyle {\mathbf{ X}}} $に変換
- 構文解析に失敗 (MathML、ただし動作しない場合はSVGかPNGで代替(最新ブラウザーや補助ツールに推奨): サーバー「https://en.wikipedia.org/api/rest_v1/」から無効な応答 ("Math extension cannot connect to Restbase."):): {\displaystyle \mathbf{Y} = \boldsymbol{\mu} + A \mathbf{X} }
を施せばよい. ただし, ベクトルはすべて列ベクトルとし, $構文解析に失敗 (MathML、ただし動作しない場合はSVGかPNGで代替(最新ブラウザーや補助ツールに推奨): サーバー「https://en.wikipedia.org/api/rest_v1/」から無効な応答 ("Math extension cannot connect to Restbase."):): {\displaystyle A}
$は$$($構文解析に失敗 (MathML、ただし動作しない場合はSVGかPNGで代替(最新ブラウザーや補助ツールに推奨): サーバー「https://en.wikipedia.org/api/rest_v1/」から無効な応答 ("Math extension cannot connect to Restbase."):): {\displaystyle A^\top}
$は$構文解析に失敗 (MathML、ただし動作しない場合はSVGかPNGで代替(最新ブラウザーや補助ツールに推奨): サーバー「https://en.wikipedia.org/api/rest_v1/」から無効な応答 ("Math extension cannot connect to Restbase."):): {\displaystyle A}
$の転置行列)を満たす下三角行列であり, コレスキー(Cholesky)分解によって求める.
[アーラン分布]
フェイズが$構文解析に失敗 (MathML、ただし動作しない場合はSVGかPNGで代替(最新ブラウザーや補助ツールに推奨): サーバー「https://en.wikipedia.org/api/rest_v1/」から無効な応答 ("Math extension cannot connect to Restbase."):): {\displaystyle k} $(正整数)のアーラン分布の密度関数は
- 構文解析に失敗 (MathML、ただし動作しない場合はSVGかPNGで代替(最新ブラウザーや補助ツールに推奨): サーバー「https://en.wikipedia.org/api/rest_v1/」から無効な応答 ("Math extension cannot connect to Restbase."):): {\displaystyle f(x)=\lambda \mathrm{e}^{-\lambda x}(\lambda x)^{k-1}/(k-1)! \ \ \ \ \ \ \ (x\geq 0) }
である. これは, パラメータ(平均値の逆数)が$構文解析に失敗 (MathML、ただし動作しない場合はSVGかPNGで代替(最新ブラウザーや補助ツールに推奨): サーバー「https://en.wikipedia.org/api/rest_v1/」から無効な応答 ("Math extension cannot connect to Restbase."):): {\displaystyle \lambda}
$の指数分布に従う$構文解析に失敗 (MathML、ただし動作しない場合はSVGかPNGで代替(最新ブラウザーや補助ツールに推奨): サーバー「https://en.wikipedia.org/api/rest_v1/」から無効な応答 ("Math extension cannot connect to Restbase."):): {\displaystyle k}
$個の確率変数の和の分布であるから,
とすればよい.
[ポアソン分布]
平均値が$構文解析に失敗 (MathML、ただし動作しない場合はSVGかPNGで代替(最新ブラウザーや補助ツールに推奨): サーバー「https://en.wikipedia.org/api/rest_v1/」から無効な応答 ("Math extension cannot connect to Restbase."):): {\displaystyle \lambda} $のポアソン乱数は, 平均が$構文解析に失敗 (MathML、ただし動作しない場合はSVGかPNGで代替(最新ブラウザーや補助ツールに推奨): サーバー「https://en.wikipedia.org/api/rest_v1/」から無効な応答 ("Math extension cannot connect to Restbase."):): {\displaystyle 1/\lambda} $の指数乱数との関係を利用して作ることができる. すなわち, $構文解析に失敗 (MathML、ただし動作しない場合はSVGかPNGで代替(最新ブラウザーや補助ツールに推奨): サーバー「https://en.wikipedia.org/api/rest_v1/」から無効な応答 ("Math extension cannot connect to Restbase."):): {\displaystyle U_1U_2\cdots U_n>{\mbox{\rm e}}^{-\lambda}} $が成り立つ最大の整数$構文解析に失敗 (MathML、ただし動作しない場合はSVGかPNGで代替(最新ブラウザーや補助ツールに推奨): サーバー「https://en.wikipedia.org/api/rest_v1/」から無効な応答 ("Math extension cannot connect to Restbase."):): {\displaystyle n} $を$構文解析に失敗 (MathML、ただし動作しない場合はSVGかPNGで代替(最新ブラウザーや補助ツールに推奨): サーバー「https://en.wikipedia.org/api/rest_v1/」から無効な応答 ("Math extension cannot connect to Restbase."):): {\displaystyle X} $とすればよい. 1個のポアソン乱数を発生するのに必要な一様乱数の個数の平均値は$構文解析に失敗 (MathML、ただし動作しない場合はSVGかPNGで代替(最新ブラウザーや補助ツールに推奨): サーバー「https://en.wikipedia.org/api/rest_v1/」から無効な応答 ("Math extension cannot connect to Restbase."):): {\displaystyle \lambda+1} $であるから, $構文解析に失敗 (MathML、ただし動作しない場合はSVGかPNGで代替(最新ブラウザーや補助ツールに推奨): サーバー「https://en.wikipedia.org/api/rest_v1/」から無効な応答 ("Math extension cannot connect to Restbase."):): {\displaystyle \lambda} $が大変に大きいときには, この方法は効率が悪いので, 二者択一法を使うほうがよい.
参考文献
[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.