乱数生成
はじめに
擬似乱数が生成できるということは，事象のシミュレーション，確率およびその他の量の推定，ランダムな割当てや選択，記号的結果の数値的な検証等にとって重要である．このような応用分野には，一様に分布した数，非一様に分布した数，復元抽出法による元，非復元抽出法による元が必要となる場合がある．
関数
RandomReal，
RandomInteger，
RandomComplexは，一様分布の乱数を生成する．
RandomRealと
RandomIntegerは組込み分布の数も生成する．
RandomPrimeはある範囲内の素数を生成する．関数
RandomChoiceと
RandomSampleはリストから復元抽出あるいは非復元抽出を行う．元は等価の重みを持つことも，非等価の重みを持つこともある．乱数生成の付加的なメソッドおよび分布を定義するためのフレームワークも含まれている．
一連の再現しない事象は，
RandomSampleでシミュレートすることができる．例えば，1から
n までの整数を順に無作為抽出する確率をシミュレートする．
これは， n が2から8ののときに， n 個の元を順に得る確率を推定する．
| Out[1]= |  |
|
| Out[2]= |  |
|
乱数生成は，モンテカルロ法の中核にある．関数
f の期待値の推定は，所望の分布からの値を生成し，その値に適用される
f の平均を見付けることにより得られる．
| Out[3]= |  |
|
この場合，すいては厳密な結果と比較することができる．
| Out[4]= |  |
|
乱数処理は，所望の特性を持つ一連の数を生成することによりシミュレートすることができる．ランダムウォークは，擬似乱数を再帰的に足すことにより生成することができる．
| Out[5]= |  |
|
乱数の代入は，記号的な式の等価性を検証するために使用することができる．例えば，2つの式の差分絶対値はランダムに生成された点で評価し，式の不等性を検証することができる．
次は，実数値のとき  と  x は異なるという証拠にはならない．
| Out[6]= |  |
|
次は，少なくともいくつかの複素数値のときは，  と xが異なるという証拠となる．
| Out[7]= |  |
|
RandomPrimeは等確率で素数を選ぶ．これは，例えばRSA暗号のために大きい素数を生成するのに便利である．素数は範囲内の素数上には一様に分布するが，全体の範囲では一様に分布しない．これは素数は一般に正の整数の範囲上で一様に分布しないためである．
| Out[8]= |  |
|
乱数生成関数
乱数
RandomRealは指定された実数値範囲上で実数の擬似乱数を生成する．
RandomIntegerは指定された整数値範囲上で整数の擬似乱数を生成する．
RandomComplexは複素平面の指定された矩形領域上で複素数の擬似乱数を生成する．
RandomPrimeは範囲内で等しい確率で素数を生成する．
実数乱数の生成
整数乱数の生成
複素数乱数の生成
素数乱数の生成
領域が
xmin および
xmax で指定されている場合，
RandomRealおよび
RandomIntegerは指定された範囲に一様に分布する数を生成する．領域が分布で指定されている場合は，その分布に対して定義される規則が使われる．また，
新しいメソッドと分布を定義するためのメカニズムも含まれている．
2引数インターフェースは，複数の乱数をすぐに得るのに便利である．さらに重要なことは，多数の擬似乱数を一度に非常に効率よく生成できるという利点がある．
0と1の間の 107個の数の生成には数分の1秒しかかからない．
| Out[9]= |  |
|
1回に1個ずつ 107個の数を生成すると，ほぼ5倍の時間がかかる．
| Out[10]= |  |
|
次元
n1 から
nk までの多次元配列では，要求される擬似乱数の合計
ntotal=
ni が生成されてから分割される．これによりできるだけ効率的に多次元配列が生成できる．それは，乱数値の合計数は可能な限り効率的に生成され，分割に必要な時間はほんのわずかしか必要ないからである．
100×100×100×10配列に必要な時間は， 107個の数のベクトルに必要な時間にほぼ等しい．
| Out[11]= |  |
| Out[12]= |  |
|
同じ次元の配列を1度に10個の数で生成すると，数倍の時間がかかる．
| Out[13]= |  |
|
統計分布では，1度に多くの数を生成することによる速度的な利点がさらに大きい．使用される同じ数の生成器から継承される効率面での利点に加え，統計分布の多くには，初等関数および特殊関数の評価のベクトル化からの利点もある．例えば，
WeibullDistributionでは初等関数
Power，
Times，
Logのベクトル評価が有利に働く．
105個のワイブル数の生成には，実質的に時間がかからない．
| Out[14]= |  |
|
105個のワイブル数を1度に1つずつ生成すると，数秒かかる．
| Out[15]= |  |
|
乱数生成は実地調査で便利なことがある．例えば，長い桁の数列の中で，ランダムな数列を探す場合等である．
これにより5つのランダムな10進数字を文字列に変換する．
| Out[16]= |  |
|
次は  の最初の100万桁を整数の文字列にする． |
これで， の最初の100万桁で先の10進数字が出現する場所が分かる．
| Out[18]= |  |
|
乱数生成は，閉形式の結果が既知でない，あるいは計算的に難しいことが分かっている分布の推定におい手も非常に便利である．乱数行列の特性でその一例を示す．
下は，一様分布の実数の5×5行列が実数の固有値を持つ確率を推定する．
| Out[19]= |  |
|
次は標準の正規分布の数の行列について同じことを行う．
| Out[20]= |  |
|
多変量分布をシミュレートする例は，ベイズ統計で使用されるギブスサンプラーである[
1]．ギブスサンプラーは，他の座標を条件とする各座標の分布が既知であるならば，多変量分布からの値をシミュレートする方法を与える．ある制約条件下では，条件付き分布から繰返し抽出することにより構築された乱数ベクトルの分布が，本当の多変量分布に収束することがある．
次の例はCasellaとGeorgeの例[
2]に対するギブスサンプラーを構築する．関心のある分布は2変量である．
y が与えられたときの
x の条件付き分布は二項分布，
x が与えられたときの
y の条件付き分布はベータ分布である．CasellaとGeorgeが述べているように，ギブスサンプラーを使った収束の検出および抽出のさまざまな方法が提案されている．簡単に言うと，収束は反復1000回以内に起ると仮定する．分布からの大きさ
n の標本は，1000回目の反復の後，
n 値とされる．しかし，この
n 値は従属値である．
二項およびベータの条件付き分布でサンプラーを定義する． |
ギブスサンプラーは，乱数生成の分布フレームワーク内部の分布オブジェクトとして定義することもできる．この特定の
分布オブジェクトとしてのギブスサンプラーの例は，
「分布の定義」に記述してある．
| Out[25]= |  |
|
第2座標の周辺分布は，ヒストグラムで可視化できる．
| Out[27]= |  |
|
条件付き分布は，それに対する十分なデータがある場合，仮定される二項およびベータ分布に近くマッチしなければならない．周辺分布の密度が高いときに最大量のデータが生じるので，それらの値は比較に使うことができる．次のグラフは，連続値の確率推定に幅.05のビンを使って，実際の条件付き分布と仮定された条件付き分布とを比較する．
0.3`≤y<0.35`の範囲での x の実際の分布と理論的分布を比較する．
| Out[28]= |  |
|
x=1のときの y の実際の分布と理論的分布を比較する．
| Out[29]= |  |
|
任意精度の実数と複素数
RandomRealおよび
RandomComplexは，デフォルトで機械精度数を生成する．任意精度数は
WorkingPrecisionオプションを設定して得ることができる．
RandomRealおよびRandomComplexのオプション
このオプションは一様分布の実数および組込み分布からの複素数・実数に有効である．
WorkingPrecisionはユーザ定義の分布に組み込むこともできる．
| Out[30]= |  |
|
| Out[31]= |  |
|
WorkingPrecisionを増加すると，精度の損失が予想されたり非常に正確な結果が必要とされるシミュレーションに便利である．精度を増加すると，計算における精度の損失を推定することもできる．
区間 [0, 1000]での J1の計算における最悪の精度の損失を推定する．
| Out[32]= |  |
|
入力の精度が指定された
WorkingPrecisionよりも小さい場合，関数により問題が警告される．入力の精度は，所望の精度の擬似乱数を生成するために人工的に増加される．
機械数7.5の精度は50より小さいので，警告が表示される．
| Out[33]= |  |
|
WorkingPrecisionは
RandomIntegerのオプションではない．整数は無限精度を持つので，精度は関数名で完全に指定される．
| Out[34]= |  |
|
ランダムな要素
RandomChoiceと
RandomSampleは，可能な要素のリストから擬似乱数的に選択する．要素は数値でも非数値でもよい．
| RandomChoice[{e1,e2,...}] | ei のひとつを擬似乱数的に選択する |
| RandomChoice[list,n] | list からn 個の要素を擬似乱数的に選択したリストを返す |
| RandomChoice[list,{n1,n2,...}] | list から擬似乱数的に選択したn1×n2×...を返す |
| RandomChoice[{w1,w2,...}->{e1,e2,...}] |
| wi により重みの付いた擬似乱数的に選択した要素を返す |
| RandomChoice[wlist->elist,n] | n 個の重み付きの擬似乱数的に選択した要素のリストを返す |
| RandomChoice[wlist->elist,{n1,n2,...}] |
| 重み付きの選択要素のn1×n2×...配列を返す |
リストからのランダムな選択
| RandomSample[{e1,e2,...},n] | ei の中からn 個の無作為抽出を行う |
| RandomSample[{w1,w2,...}->{e1,e2,...},n] |
| 重みwiを使って選ばれたei の中からn 個の無作為抽出を行う |
| RandomSample[{e1,e2,...}] | ei の擬似ランダム置換を行う |
| RandomSample[wlist->elist] | 最初の重みwlist を使ってelist の擬似ランダム置換を行う |
リストからの無作為抽出
RandomChoiceと
RandomSampleの主な違いは，
RandomChoiceが
ei から復元抽出するのに対し，
RandomSampleは非復元抽出するという点である．
RandomChoiceで選ばれる要素の数は
elist の要素の数で制限されず，要素
ei は複数回選択されてもよい．
RandomSampleで返される標本の数は
elist の要素の数で制限され，その標本の中で異なる要素が現れる回数は，
elist でその要素が現れる回数で制限される．
RandomChoiceまたは
RandomSampleの第1引数がリストならば，要素は同じ確率で選ばれる．重み指定は，
ei の集合の分布を定義する．この重みは正でなければならないが，和が1である必要はない．重み
{w1, ..., wn}では，初期分布の
ei の確率は
wi/
wjである．
RandomSampleは非復元抽出を行うので，各選択の後，残りの合計の重みに基づいて重みが内部的に更新される．
RandomChoiceは可能な結果の有限リストを持つ，独立で同一の分布の事象のシミュレーションに使うことができる．
| Out[35]= |  |
|
| Out[36]= |  |
|
RandomChoiceは有限サポートを持つどの離散分布から観測を生成するためにも使うことができる．
| Out[37]= |  |
|
シミュレートされた1000個の点に対する経験的確率密度関数である．
| Out[40]= |  |
|
RandomSampleは，結果のリストの各要素が1度だけしか観測されないような結果の有限集合からの観測をシミュレーションするために使うことができる．そのリストには異なる値が複数回出現してもよい．
次は青が80個，赤が45個入った容器から7回引くことをシミュレートする．
| Out[41]= |  |
|
リストのすべての要素を無作為に抽出すると，ランダム置換になる．
| Out[42]= |  |
|
要素に重みを割り当てると，置換において大きい重みを持つ値が，小さい値を持つ重みよりも先に現れる傾向にあるランダム置換となる．
| Out[43]= |  |
|
重み付きあるいは重みなしの要素を持つ同じリストでは，
RandomSample[#, 1]&は分布が
RandomChoiceに等しい．
大きさ1の 105個の無作為抽出要素に対する経験的確率密度関数．
| Out[45]= |  |
| Out[46]= |  |
|
| Out[47]= |  |
|
この2つの例に対する確率は互いに，また理論的値に非常に近い．
| Out[48]= |  |
|
RandomSampleは，臨床試験等における，群への無作為割付けに使うこともできる．次では整数を使っているが，名前やID番号等の識別値を使うこともできる．
下は，20の元を同じ大きさの4つの群に無作為に割り付ける．
| Out[49]= |  |
|
RandomChoiceおよび
RandomSampleは
SeedRandomの
Methodオプションの変更の影響を受けることがある．組込みメソッドは
「メソッド」で説明する．また，新しいメソッドを定義するメカニズムは
「自分の生成器を定義する」で説明する．
種と局所化
擬似乱数生成では，明らかなランダム性のレベルを持つ数がアルゴリズムにより生成される．擬似乱数生成のためのメソッドでは，通常循環関係を使って現在の状態から数を生成し，次の数を生成する新しい状態を構築する．その状態は，アルゴリズムの循環関係を初期化するために使われる整数を生成器に撒くことで設定することができる．
「種」と呼ばれる初期値が与えられると，擬似乱数生成器は完全に確定的になる．多くの場合，乱数生成に対する種を局所的あるいは大域的に設定し，「乱数」値の数列を継続的に得ることが望ましい．大域的に設定すると，新しい種を明示的に設定しない限り，この種が今後の擬似乱数に影響する．局所的に設定すると，種は局所化されたコード内の乱数と要素の生成にのみ影響を与える．
局所化と種関数
SeedRandom関数は乱数生成器に種を蒔く方法を提供する．
SeedRandomだけを使うと，乱数生成器用の種は大域的に設定される．
BlockRandom関数は，大域的な状態に影響を与えずに乱数生成器用の種を局所的に設定あるいは変更する方法を提供する．
| Out[50]= |  |
|
次では最初の
RandomRealが
BlockRandom内で生成され，2つ目が
BlockRandomの外側で生成されているため，2つの異なる数を返す．
| Out[51]= |  |
|
SeedRandomは乱数生成を切り替えるメカニズムも提供する．
SeedRandomのオプション
個々の生成器には
Methodオプションでその生成器を指定することにより，直接種を蒔くことができる．すべての生成器には，設定
Method->Allで種を蒔くことができる．
ここではデフォルトの生成器の種は 1であるが， "Rule30CA"生成器には種がない．
| Out[52]= |  |
|
"Rule30CA"に種 1を蒔くと，異なる乱数が生成される．
| Out[53]= |  |
|
メソッド
すべてのシステムで，5つの擬似乱数生成メソッドが使用できる．6番目のプラットフォーム依存のメソッドはIntelベースのシステムで利用できる．
「自分の生成器を定義する」で説明してある新しいメソッド定義のためのフレームワークも含まれている．
| "Congruential" | 線形合同法（低質のランダム性） |
| "ExtendedCA" | 拡張セルオートマトン法（デフォルト） |
| "Legacy" | Mathematica 6.0より前のデフォルト法 |
| "MersenneTwister" | メルセンヌ・ツイスタシフトレジスタ法 |
| "MKL" | インテルMKL（インテルベースのシステム） |
| "Rule30CA" | ウルフラム則30生成法 |
組込みメソッド
シードが2020のときの各メソッドからの整数の擬似乱数．
| Out[54]= |  |
|
| Out[55]= |  |
|
Congruential
"Congruential"は線形合同法を使う．これは擬似乱数が
xi/m から得られる0と1である，最も単純な擬似乱数生成法のひとつである．
xi は下の漸化式
であり，固定された整数
b，
c，
m はそれぞれ乗数，増分，法と呼ばれる．増分が0ならば，生成器は乗法合同法である．
b，
c，
m の値はオプションを使って
"Congruential"メソッドに設定することができる．
| | |
| "Bits" | Automatic | ビットから構成される数に使用するビットの範囲を指定する |
| "Multiplier" | 1283839219676404755 | 乗数の値 |
| "Increment" | 0 | 増分の値 |
| "Modulus" | 2305843009213693951 | 法の値 |
| "ConvertToRealsDirectly" | True | 合同関係から直接実数を構築するかどうか |
Method "Congruential"のオプション
線形合同法は周期的であり，多数の乱数値が必要な場合は特にランダム性の質が低下する傾向にある．実数が合同関係から直接生成される場合，周期は
m以下である．
デフォルトのオプション値は，大きい周期を持ち，64ビット効率に対応するよう選ばれる．デフォルトオプションでは，
"Congruential"生成器は，合同数生成に関する継承問題があるにもかかわらず，乱数性をみる多数の標準テストに合格する．
乗法合同法の周期は，法に互いに素である法と等しいかそれ以下の正の整数の数により制限される．この上限は，オイラーのトーティエント関数である．
| Out[57]= |  |
|
実際の周期は
i
bi mod
m となるような最小の整数
i を見付けることで決めることができる．
| Out[58]= |  |
|
データを6つの要素の集合に分割すると再帰が分かる．
| Out[59]= |  |
|
生成された数列をプロットすると，異なる数字がグラフィカルに表示される．
| Out[60]= |  |
|
"ConvertToRealsDirectly"が
Falseに設定されると，数列の要素から1度に8ビット取り52ビットの機械精度数を構築することにより実数が生成される．この方法で生成された合同数はまだ循環するが，循環は初期の合同関係の反復ではなくビットパターンの反復に依存する．
"Bits"オプションは
Automatic，非零の整数，ビットからの数を構築するために使われる法
m のビット範囲を指定する2つの非零の整数のリストのいずれかを取る．
Automaticでは
{2, -1}が使われるが，
m が2のベキ乗のときは
{1, -1}が使われる．
ExtendedCA
デフォルトの
"ExtendedCA"メソッドは，セルオートマトンを利用して高質の擬似乱数を生成する．この方法は，特殊な5近傍則を使うので，それぞれの新規セルは，直前のステップの5つの非隣接セルに依存する．
セルオートマトンベースの乱数生成は，決定的規則に従って0と1の状態ベクトルを進化させる．指定されたセルオートマトンでは，新しい状態ベクトルの指定の位置にある要素（またはセル）は，古い状態ベクトルにおけるそのセルの一定の隣接セルにより決められる．状態ベクトルのセルの部分集合は，ランダムなビットとして出力され，それから擬似乱数が生成される．
"ExtendedCA"で使用されるセルオートマトンは，非常に高レベルのランダム性を生成する．そのランダム性は非常に高度なので，出力ですべての単独セルを使った場合でさえも，1つのセルと5つの前のセルの間に明らかな相関があるにもかかわらず，多数のランダム性テストにパスするようなビットのストリームを生成する．
状態ベクトルの大きさおよび飛ばされるセルを変更するための2つのオプションが含まれている．デフォルトは質と速さで選ばれ，通常このオプションを変更する必要はない．
| | |
| "Size" | 80 | 64の乗数としての状態ベクトルの大きさ |
| "Skip" | 4 | 飛ばすセルの数 |
Method "ExtendedCA"のオプション
使用される状態ベクトルの長さはデフォルトで80×64=5120セルに設定されている．64の乗数が
"Size"オプションで指定できる．
実際には，各状態ベクトルで4つおきにセルを使うことは，非常に厳しいランダム性テストにパスするのに十分であると分かっている．これは
"Skip"オプションに使われるデフォルトである．さらに高速な乱数生成法では，
"Skip"設定が2もしくは1にもなるが，乱数の質は低下する．
"ExtendedCA"はデフォルトの数生成法である．
| Out[61]= |  |
| Out[62]= |  |
|
Legacy
"Legacy"メソッドは，
Mathematica バージョン6.0より前では
Randomにより呼び出される生成法を使用する．実数にはマーサグリア・ザマンボロー付き減算ジェネレータが使われる．整数の生成法にはウルフラム則30のセルオートマトンジェネレータを使う．規則30ジェネレータは，小さい整数の場合は直接，大きい整数の場合は一定のビットを生成するために使われる．
| Out[63]= |  |
|
| Out[64]= |  |
|
バージョン6.0以前に生成された列との一貫性を保障するため，
Automaticメソッドの種の集合が
"Legacy"メソッドにも適用される．
MersenneTwister
"MersenneTwister"はMatsumotoとNishimura [
3][
4]によるメルセンヌツイスタ生成器を使用する．メルセンヌツイスタは，周期
219937-1の一般化されたフィードバックシフトレジスタ生成器である．
| Out[65]= |  |
|
"MersenneTwister"メソッドにはオプションがない．
MKL
"MKL"メソッドは，IntelのMKLライブラリで提供されている乱数生成器を使う．MKLライブラリは，プラットフォーム依存である．
"MKL"メソッドは，Intelベースのシステムでのみ使用可能である．
Method "MKL"のオプション
| "MCG31" | 31ビット乗算合同法 |
| "MCG59" | 59ビット乗算合同法 |
| "MRG32K3A" | 3次の2つの要素を持つ結合multiple recursive generator |
| "MersenneTwister" | メルセンヌツイスタシフトレジスタ生成器 |
| "R250" | 一般化フィードバックシフトレジスタ生成器 |
| "WH" | Wichmann-Hillの結合乗算合同法 |
| "Niederreiter" | Niederreiterの低食い違い量列 |
| "Sobol" | Sobolの低食い違い量列 |
"MKL"メソッド
最初の6つの方法は一様乱数生成法である．
"Niederreiter"と
"Sobol"はNiederreiter列と Sobol列を生成する．これらの数列は非一様であり，数値法で便利な場合もある内在構造を持つ．例えば，これらの数列は通常多次元のモンテカルロ積分で速く収束する．
| Out[66]= |  |
|
| Out[67]= |  |
|
Rule30CA
"Rule30CA"メソッドはウルフラム則30のセルオートマトン生成器を使う．以下の関係
を使い，0と1の状態ベクトルを進化させることによりビットを得る．ここで
f (i, t)は時間
tにおけるセル
i の値である．
| | |
| "Size" | 9 | 29の乗数としての状態ベクトルの大きさ |
Method "Rule30CA"のオプション
使用される状態ベクトルの長さは，デフォルトで
9×29=261セルに設定されている．29の乗数は
"Size"オプションで指定することができる．
"Rule30CA"を使って整数の乱数の2x3x4テンソルを与える．
| Out[68]= |  |
|
"Rule30CA"メソッドは，各状態ベクトルから最初のビットしか使わないので，各状態ベクトルから複数のビットを使用する
"ExtendedCA"メソッドよりも遅い．
自分の生成器を定義する
メソッドは正しいテンプレートに従っている限り，乱数フレームワークにプラグインすることができる．生成期のオブジェクトは
gsym[data]という形式である．ここで
gsym は生成器を同定し，規則が接続されているシンボルである．
data は実質的に，生成器の定義に関連付けられたトップレベルの評価に対してプライベートである．
生成器の初期化は
Random`InitializeGeneratorの呼出しにより実行される．
| Random`InitializeGenerator[gsym,opts] |
| 生成器gsym をオプションoptsで初期化する |
生成器初期化関数
Random`InitializeGeneratorは形式
gsym[data]の生成器オブジェクト
gobj を返すことが想定されている．
生成器はランダムなビットストリーム，ランダムな整数，ランダムな実数の生成をサポートすることができる．生成器がビットストリームをサポートする場合，実数と整数はビットストリームの変換により生成できる．メソッドの設定時には，何をどのようにサポートするかを決定するための属性が問い合される．
| GeneratesBitsQ | メソッドがビットを生成する場合はTrueに設定 |
| GeneratesIntegersQ | メソッドが指定された範囲の整数を生成するならばTrueに設定 |
| GeneratesRealsQ | メソッドが指定された範囲と精度で実数を生成するならばTrueに設定 |
生成器の属性
ビットストリームがサポートされる場合，
gobj["GenerateBits"[nbits]]は
n 個のランダムなビット，あるいは項目が0か1の長さ
nbits のリストを返すことが想定される．
ランダムな整数がサポートされる場合，
gobj["GenerateIntegers"[n, {a, b}]]は範囲
[a, b]の
n 個のランダムな整数のリストを返すことが想定される．結果が範囲を超えると警告メッセージが出される．
ランダムな実数がサポートされる場合，
gobj["GenerateReals"[a, {a, b}, prec]]は範囲
[a, b]で精度
prec の
n 個の実数のリストを返すことが想定される．結果が範囲を超えたり，異なる精度であったりした場合は，警告メッセージが出される．
生成関数のすべてに対して，返されるものは
{res, gobj}である．
res は正しいタイプの結果であり，
gobj は新しい生成器オブジェクト（状態の変化を反映して）である．
整数
seed に対する種まきは
gobj["SeedGenerator"[seed]]により実行される．
gobj["SeedGenerator"[seed]]は新しい生成器オブジェクトを返すものと想定される．
例題：乗算合同法
次の例では，乗算合同法を定義する．乗算合同法は，漸化式
に従う．以下で定義されるように，生成器は実数の生成だけができる．
生成器 MultiplicativeCongruentialのデフォルトオプションを設定する． |
生成器の初期化は乗数と法の値を抽出する．この値のいずれかが正の整数でない場合，初期化は失敗する．
MultiplicativeCongruentialが実数を生成するよう設定する． |
実数生成器は所望の数の実数
n 個と新しい
MultiplicativeCongruential生成器を返す．新しい生成器の種は，漸化式に基づき更新される．
MultiplicativeCongruential生成器を使い，10個の実数を生成する．
| Out[74]= |  |
|
| Out[75]= |  |
|
例題：Blum-Blum-Shub法
Blum Blum Shub法は暗号化目的の擬似乱数ビットを生成するための2次合同法である[
5]．合同式は指定された素数
p と
q に対する法
p×
q である．
生成器 BlumBlumShubのデフォルトオプションを指定する． |
以下で補助関数と生成器に対するエラーメッセージを定義する． |
生成器の初期化でオプション値が抽出され，実際の生成器を呼び出す前に，必要に応じてエラーメッセージが出される．
BlumBlumShubがビット生成器であることを指定し，ビット幅を決める． |
BlumBlumShub生成器を使い，5つの整数と5つの実数を生成する．
| Out[86]= |  |
|
統計分布
非一様分布からランダムな変量を生成する目的は，0と1の間のランダムな一様変量を生成し，所望の分布の乱数値の累積分布関数の逆関数を計算することにある．しかし，実際には，直接これに従うと多数の乱数変量が必要な場合に計算的に非常に厳しくなる．特に累積分布関数の逆関数が複雑で，閉形式で表せない場合はそうである．
そのような場合，累積分布関数の直接の逆関数にはテーブル検索，分布関係に基づく直接構造，棄却採択法の方が効率的であることが多い．あるレベルにおいては，このような方法はすべて一様分布に従う
RandomReal値，非一様分布の
RandomInteger値，重み付きの
RandomChoiceの観測値，あるいはこれらすべての値の組合せに依存する．その結果，
SeedRandomで設定された方法は，統計分布の乱数観測値に影響を及ぼす．
Mathematica の多数の分布に対して
RandomRealおよび
RandomIntegerで使用される方法はGentle [
6]で提唱されるか記述されている方法に従っており，バージョン6.0より前の
Mathematica に含まれていた標準アドオンの
Randomおよび
RandomArrayで使用された方法と必ずしも同じではない．
すべての組込みの統計分布からの乱数観測値は，
RandomRealか
RandomIntegerを使って生成することができる．
連続分布からの乱数値の生成
離散分布からの乱数値の生成
1変量および多変量の連続分布からの観測は
RandomRealで得られ，その離散分布からの観測は
RandomIntegerで得られる．もし
RandomIntegerが連続分布に使われたり
RandomRealが離散分布に使われたりすると，エラーが出され，入力は評価しないままで返される．
 2分布は連続分布なので，次は失敗する．
| Out[87]= |  |
|
WorkingPrecisionは，範囲上の一様の数に対するのと同様に，連続分布に対する
RandomRealのオプションである．
| Out[88]= |  |
|
多変量分布および多変量正規分布から派生する分布は多変量統計パッケージに含まれている．これらの分布に対する生成器は，パッケージをロードしてアクセスすることができる．
| Out[90]= |  |
|
| Out[91]= |  |
|
連続分布
複数の一様変数から，あるいは一様分布以外の変数からの単独の確率変数の直接構築も使用される．正規変数は，Box-Müller法を使い一様乱数のペアから，ペアとして生成される．
HalfNormalDistributionおよび
LogNormalDistributionの変数は正規変数の直接変換により得られる．多変量統計パッケージの
MultinormalDistributionおよび
QuadraticFormDistributionも正規変数からの直接構築を使う．
InverseGaussianDistributionは，正規変数と一様変数を含むacceptance-complement法を使用する．この方法はMichael，Schucany，Haasによるものであり，Gentle [
6]に記述されている．
MaxwellDistribution変数は
ChiDistribution変数から構築される．カイ変数自体は，
GammaDistributionの特殊形である
ChiSquareDistributionから得られる．
ほとんどの場合，
FRatioDistributionは単独のベータ変数乱数からそれぞれの乱数値を構築する．自由度が小さい場合，
FRatioDistribution変数は，ベータ構築で発生する可能性のある0による割り算の可能性を避けるために，ガンマ変数のペアから生成される．
NoncentralChiSquareDistribution[
, 
]，
，の変数生成では，
2分布の加法的性質を使い，非整数
に対する累積密度関数の高価な逆関数の計算を避ける．この加法的性質は，例えば Johnson，Kotz，Blakrishnan [
7]で与えられる．
=1のとき，非心
2変数は平均
，分散1の正規変数の平方として生成することができる．
≠1では，非心
2変数は中心および非心
2確率変数の和として得られる．
>1のとき，
X+Y はX~
およびY~
ならば
分布に従う．この関係は
<1のときは使えない．その場合，構築は
および
の
X+Yとなる．ここで
は
が0に近付くにつれ極限非心
2分布となる．極限分布
はポワソン変数と
2変数を合わせたものであり，0において非零の確率質量を，正の値の時は連続密度を持つ．
NoncentralFRatioDistribution変数は1つの中心
2変数と1つの非心
2変数から得られる．
多変量統計パッケージの
WishartDistributionでは，行列はSmithとHockingの方法[
8]で生成される．この方法は，カイ分布に従った対角要素と正規分布に従った非対角要素を持つ行列からウィッシャート行列を構築する．
GammaDistribution[
, 
]指数変数は
=1のときに生成される．それ以外のときは，Cheng，Feast [
9]によるメソッドおよびAhrens，Dieter [
10]によるメソッドが使われる．
ベータ変数はベータパラメータ
および
の値により複数の方法を切り替えることにより構築される．両方のパラメータが1ならば，一様確率変数が生成される．ベータパラメータのうちのひとつが1ならば，累積密度関数の閉形式の逆関数評価が使われる．それ以外のとき，
RandomRealはJöhnk [
11]，Cheng [
12]，Atkinson [
13]による棄却採択法のいずれかを使う．2つのガンマからの構築よりも棄却採択法を使うことの利点の一例を下に挙げる．直接棄却採択法はガンマペア構築法のほぼ2倍の速さである．
| Out[92]= |  |
| Out[93]= |  |
|
StudentTDistributionでは，
RandomRealで使用されるメソッドはBailey [
14]のPolar Rejection法である．この方法は以下で示す通り，正規変数および
2 変数からの直接構築よりも効率的である．直接構築は100万のスチューデントの
t 変数に対してPolar法のほぼ1.5倍の時間がかかる．
スチューデントの t に対する直接構築とBaileyのPolar Rejection法を比較する．
| Out[94]= |  |
| Out[95]= |  |
|
離散分布
整数の乱数生成で表検索が使われるときは，重みに対する分布の確率質量関数を使い，
RandomChoiceにより実装される．ほとんどの離散分布は，所望の標本サイズが与えられたときに，重みのリストの構築が効果になると見込まれると，別の方法に切り替える．例えば，
p が1に近付くと，
LogSeriesDistribution[p]は直接構築
に切り替える．ここで，
U および
V は区間
[0, 1] [
15]上で一様分布に従う．また，パラメータと標本の大きさによっては，
NegativeBinomialDistribution[n, p]はポワソンとガンマを合わせた構築法に切り替える．これは平均がガンマ分布に従うポワソン変数[
6]である．
BinomialDistribution，
HypergeometricDistribution，
PoissonDistributionはすべて，確率分布関数値の計算のオーバーヘッドが所望の乱数値の数と比較して小さければ，密度関数からの直接抽出に依存する．確率分布関数値の直接計算においてオーバーフローやアンダーフローが起きると，棄却採択法も変数の生成が行えるので，乱数が生成できるパラメータ値の範囲が拡張される．
二項分布および超幾何分布はKachitvichyanukulおよびSchmeiserの棄却採択法に少し修正を加えたものに切り替える．BTPE (Binomial, Triangle, Parallelogram, Exponential) アルゴリズム[
16]の棄却採択法の部分に基づく二項法は，迅速な採択テストのために3つの領域を持つ区分的なMajorizing関数と，Minorizing三角関数を効率的に使う．Majorizing関数とMinorizing関数は再スケールされた二項分布の密度関数の中心付近に2つの平行四辺形のエンベロプを生成し，Majorizing関数の裾は，スケールされた二項分布の裾上に指数関数的なエンベロプを形成する．検索表を構築するよりもBTPE法を使った方が明らかによい場合の一例として，観測はほとんど必要とされていないのに検索表が大きいという場合がある．
KachitvichyanukulとSchmeiserのH2PEアルゴリズム[
17]の棄却採択法の部分に基づく超幾何法は，スケールされた超幾何分布の密度関数周囲に3つの領域を持つMajorizing関数を使う．密度関数の中心部は，長方形領域でエンベロプに入れられ，分布の裾は指数関数を限界とする．
PoissonDistributionで使用される棄却採択法はAhrensとDieter [
18]によるものである．棄却と採択は離散正規変数を使って実行され，
が増加するにつれて
PoissonDistribution[]が
に傾くことを利用している．
ZipfDistributionの乱数値は，Devroye [
15]により記述され手いる棄却採択法で生成される．このメソッドは一様変数のペアおよび基礎演算以外では
Floorと非整数ベキしか含まない検定を使い，ジップ分布値を効率的に得る．
分布の定義
分布の定義は
Random`DistributionVectorに対する規則でサポートされる．
DistributionVectorは指定された精度の数を持つ指定された長さのベクトルを返すことが想定されている．
| Random`DistributionVector[dist,n,prec] |
| 精度prec のdist からn 個の観測値を生成するための規則を定義する |
分布から乱数生成を定義するための関数
関数から乱数値を生成するための規則は，一般に分布の頭部の
TagSetで定義される．分布自体にパラメータが含まれることもある．簡単な例として，以下で区間
[-b, -a]上の一様分布を表す
NegativeOfUniform[a, b]に対する規則を定義する．
これで
NegativeOfUniformの乱数が，他の組込みの連続分布と同様に，
RandomRealで生成できるようになった．
以下は NegativeOfUniformから機械精度および精度20の数を与える．
| Out[97]= |  |
|
行列および高次元のテンソルも
RandomRealで直接生成することができる．
RandomRealは
Random`DistributionVectorに与えられる定義を使い，所望の乱数値の合計数を生成し，それを指定された次元に分割する．
NegativeOfUniform数の3x4行列である．
| Out[98]= |  |
|
離散分布も同様に定義することができる．大きく異なる点は，
Random`DistributionVectorの精度引数が
Infinityであるということである．
NegativeOfUniformの離散型で，簡単な例を示す．
NegativeOfDiscreteUniformの乱数値は，
RandomIntegerから得られるようになった．
NegativeOfDiscreteUniformの10個の数である．
| Out[100]= |  |
|
前の例では分布の定義に対する基本的なフレームワークを示したが，分布自体は特に面白いものではない．実際，この2つの場合では新しい分布のシンボルに定義を加える代りに，
RandomRealおよび
RandomIntegerから値だけを生成し，最終結果に-1を掛けた方が簡単である．次の例では分布に定義を加えた方が便利であるような，少し複雑な分布を示す．
例題：逆関数による正規分布
一般の1変量分布から乱数値を生成する教科書的な定義には下の2つのステップが含まれる．
- 関心のある分布に対するq の累積分布関数の逆関数を計算する
このプロセスを例示するために，ここでは正規分布を使う．この方法を使ってランダムな正規変数を生成するためには，点
q において正規分布に対する
Quantileから開始し，
q を0と1の間の一様乱数で置き換えるとよい．
| Out[101]= |  |
|
これで新しい分布オブジェクトを使い，累積分布関数の逆関数を使う正規乱数生成器を定義することができる．
次は逆関数により生成された10個の正規乱数である．
| Out[103]= |  |
|
標本平均および標準偏差を伴う 104個の正規乱数の標本である．
| Out[104]= |  |
| Out[105]= |  |
|
正規分布は，直接倒置以外の方法が好まれる場合のひとつである．累積分布関数の逆関数が擬似正規乱数生成に対して完全に有効なメソッドである反面，完全に効率的であるわけではない．
InverseErfの数値評価は，
NormalDistributionによって使用されるBox-Müller法で要求される正弦関数や対数関数の評価よりも計算的に効果である．
組込みの方法では同じ数の値に対してほとんど時間がかからない．
| Out[106]= |  |
| Out[107]= |  |
|
例題：円板上の一様分布
Random`DistributionVectorは多次元分布の生成器を定義するために使うこともできる．例えば，単位円板上の一様分布からの乱数点があるとすると，
である実数点
{x, y}の集合が望まれる．そのような乱数点は次のように構築することができる．
- [0, 2)で一様分布するランダムな角度
を生成する
- [0, 1]で一様分布するランダムベクトルu を生成する
を返す
半径
r の円板上で一様分布する点を生成するために，戻った順序対に
rを掛けることができる．
| Out[110]= |  |
|
次は，この円板上の 104個の生成点の分布を可視化する．
| Out[111]= |  |
|
例題：ディリクレ分布
n 次元ディリクレ分布は正の値のベクトル
{1, ..., n}でパラメータ化され，1から
n までで
i に対して
xi=1であり
xi
[0, 1]であるようなベクトル集合
{x1, ..., xn}に対するサポートを持つ．従って，ディリクレ分布は
n 次元単位超立方体
[0, 1]nの
n-1次元の部分空間上で定義される．ディリクレ分布はベータ分布の多変量拡張である．
y が
BetaDistribution[, ]に従うならば，
{y, 1-y}はパラメータベクトル
{, }のディリクレ分布に従う．
n 次元ディリクレ変数は，
n 個のガンマ変数から生成することができる．パラメータベクトルが
{1, ..., n}のとき，プロセスは下の通りである．
以下はシンボル DirichletDistributionに付属したディリクレ生成器である． |
| Out[113]= |  |
|
例題：ギブスサンプラー
ギブスサンプラーは分布として定義することもできる．例として，ベータ分布と二項分布を組み合せたギブスサンプラーを見てみる．このサンプラーの特別な場合は
既出の例で説明してある．ここでは，分布は2つのパラメータ
m および
で定義される．
以下はギブスサンプラー BinomialBetaSamplerを定義する． |
既出の特別なギブスサンプラーでは，
m が16で
が2であった．
以下は m=16， =2のサンプラーからの5つのベクトルである．
| Out[115]= |  |
|
参考文献
[1] Geman S. and D. Geman. "Stochastic Relaxation, Gibbs Distributions, and the Bayesian Restoration of Images" IEEE Transactions on Pattern Analysis and Machine Intelligence 6, no. 6 (1984): 721-741
[2] Casella G. and E. I. George. "Explaining the Gibbs Sampler" The American Statistician 46, no. 3 (1992): 167-174
[3] Matsumoto M. and T. Nishimura. "Mersenne Twister: a 623-Dimensionally Equidistributed Uniform Pseudorandom Number Generator" ACM Transactions on Modeling and Computer Simulation 8, no. 1 (1998): 3-30
[4] Nishimura T. "Tables of 64-Bit Mersenne Twisters" ACM Transactions on Modeling and Computer Simulation 10, no. 4 (2000): 348-357
[5] Junod P. "Cryptographic Secure Pseudo-Random Bits Generation: The Blum-Blum-Shub Generator" August 1999. http://crypto.junod.info/bbs.pdf
[6] Gentle J. E. Random Number Generation and Monte Carlo Methods, 2nd ed. Springer-Verlag (2003)
[7] Johnson N. L., S. Kotz, and N. Balakrishnan. Continuous Univariate Distributions, Volume 2, 2nd ed. John Wiley & Sons (1995)
[8] Smith W. B. and R. R. Hocking. "Algorithm AS 53: Wishart Variate Generator" Applied Statistics 21, no. 3 (1972): 341-345
[9] Cheng R. C. H. and G. M. Feast. "Some Simple Gamma Variate Generators" Applied Statistics 28, no. 3 (1979): 290-295
[10] Johnson M. E. Multivariate Statistical Simulation. John Wiley & Sons (1987)
[11] Jöhnk M. D. "Erzeugung von Betaverteilten und Gammaverteilten Zufallszahlen" Metrika 8 (1964): 5-15
[12] Cheng R. C. H. "Generating Beta Variables with Nonintegral Shape Parameters" Communications of the ACM 21, no. 4 (1978): 317-322
[13] Atkinson A. C. "A Family of Switching Algorithms for the Computer Generation of Beta Random Variables" Biometrika 66, no. 1 (1979): 141-145
[14] Bailey R. W. "Polar Generation of Random Variates with the t-Distribution" Mathematics of Computation 62, no. 206 (1994): 779-781
[15] Devroye L. Non-Uniform Random Variate Generation. Springer-Verlag (1986)
[16] Kachitvichyanukul V. and B. W. Schmeiser. "Binomial Random Variate Generation" Communications of the ACM 31, no. 2 (1988): 216-223
[17] Kachitvichyanukul V. and B. W. Schmeiser. "Computer Generation of Hypergeometric Random Variates" Journal of Statistical Computation and Simulation 22, no. 2 (1985): 127-145
[18] Ahrens J. H. and U. Dieter "Computer Generation of Poisson Deviates from Modified Normal Distributions" ACM Transactions on Mathematical Software 8, no. 2 (1982): 163-179
で擬似乱数的に0か1を返す
として生成され，U はUniformDistribution[0, 1]に従う．BetaBinomialDistributionおよびBetaNegativeBinomialDistributionは，確率パラメータをBetaDistributionの確率変数とするBinomialDistribution変数およびNegativeBinomialDistribution変数から構築される．
iを生成する
