乱数生成

はじめに

擬似乱数が生成できるということは,事象のシミュレーション,確率およびその他の量の推定,ランダムな割当てや選択,記号的結果の数値的な検証等にとって重要である.このような応用分野には,一様に分布した数,非一様に分布した数,復元抽出法および非復元抽出法による元が必要となる場合がある.

関数RandomRealRandomIntegerRandomComplexは,一様分布の乱数を生成する.RandomVariateは組込み分布の数も生成する.RandomPrimeはある範囲内の素数を生成する.関数RandomChoiceRandomSampleはリストから復元抽出あるいは非復元抽出を行う.元は等価の重みを持つことも,非等価の重みを持つこともある.乱数生成の付加的なメソッドおよび分布を定義するためのフレームワークも含まれている.

一連の再現しない事象は,RandomSampleでシミュレートすることができる.例えば,1から までの整数を順に無作為抽出する確率をシミュレートする.

これは, が2から8のときに, 個の元を順に得る確率を推定する.
In[10]:=
Click for copyable input
Out[10]=
結果を,理論的な確率と比較することができる.
In[11]:=
Click for copyable input
Out[11]=

乱数生成は,モンテカルロ法の中核にある.関数 の期待値の推定は,所望の分布からの値を生成し,その値に適用される の平均を見付けることにより得られる.

次は正規分布の六次モーメントを推定する.
In[12]:=
Click for copyable input
Out[12]=
この場合,推定は厳密な結果と比較することができる.
In[13]:=
Click for copyable input
Out[13]=

乱数処理は,所望の特性を持つ一連の数を生成することによりシミュレートすることができる.ランダムウォークは,擬似乱数を再帰的に足すことにより生成することができる.

ここで0から出発するランダムウォークを生成する.
In[14]:=
Click for copyable input
Out[14]=

乱数の代入は,記号的な式の等価性を検証するために使用することができる.例えば,2つの式の差分絶対値はランダムに生成された点で評価し,式の不等性を検証することができる.

次は,実数値のときは異なるという証拠にはならない.
In[15]:=
Click for copyable input
Out[15]=
次は,少なくともいくつかの複素数値のときは,が異なるという証拠となる.
In[16]:=
Click for copyable input
Out[16]=

RandomPrimeは等確率で素数を選ぶ.これは,例えばRSA暗号のために大きい素数を生成するのに便利である.素数は範囲内の素数上には一様に分布するが,全体の範囲では一様に分布しない.これは素数は一般に正の整数の範囲上で一様に分布しないためである.

指定された範囲内の素数は等確率で生成される.
In[17]:=
Click for copyable input
Out[17]=

乱数生成関数

主な関数はRandomRealRandomIntegerRandomComplexRandomVariateRandomChoiceRandomSampleである.RandomRealRandomIntegerRandomComplexは,数値の範囲が与えられると数を発生させる.RandomChoiceRandomSampleは非数値を含むことのある有限集合から元を生成する.

乱数

RandomRealは指定された実数値範囲上で実数の擬似乱数を生成する.RandomIntegerは指定された整数値範囲上で整数の擬似乱数を生成する.RandomComplexは複素平面の指定された矩形領域上で複素数の擬似乱数を生成する.RandomVariateは指定された統計分布から擬似乱数を生成する.RandomPrimeは範囲内で等しい確率で素数を生成する.

RandomReal[]範囲0から1で実数の擬似乱数を返す
RandomReal[{xmin,xmax}]範囲 から で実数の擬似乱数を返す
RandomReal[xmax]範囲0から で実数の擬似乱数を返す
RandomReal[domain,n]n 個の実数の擬似乱数のリストを返す
RandomReal[domain,{n1,n2,}]実数の擬似乱数の ××配列を返す

実数乱数の生成

RandomInteger[{imin,imax}]範囲で整数の擬似乱数を返す
RandomInteger[imax]範囲で整数の擬似乱数を返す
RandomInteger[]確率で擬似乱数的に0か1を返す
RandomInteger[domain,n]n 個の整数の擬似乱数を返す
RandomInteger[domain,{n1,n2,}]整数の擬似乱数の ××配列を返す

整数乱数の生成

RandomComplex[]単位平方内で複素数の擬似乱数を返す
RandomComplex[{zmin,zmax}] で囲まれた長方形で複素数の擬似乱数を返す
RandomComplex[zmax]0と でかこまれた長方形で複素数の擬似乱数を返す
RandomComplex[domain,n]n 個の複素数の擬似乱数のリストを返す
RandomComplex[domain,{n1,n2,}]複素数の擬似乱数の ××配列を返す

複素数乱数の生成

RandomComplex[dist]分布 dist から擬似乱数値を返す
RandomComplex[dist,n]dist から n 個の擬似乱数値のリストを返す
RandomComplex[dist,{n1,n2,}]dist からの擬似乱数値の ××配列を返す

分布からの乱数値の生成

RandomPrime[{imin,imax}]範囲で素数の擬似乱数を返す
RandomPrime[imax]範囲2から で素数の擬似乱数を返す
RandomPrime[domain,n]n 個の素数の擬似乱数のリストを返す
RandomPrime[domain,{n1,n2,}]素数の擬似乱数の ××配列を返す

素数乱数の生成

領域が および で指定されている場合,RandomRealおよびRandomIntegerは指定された範囲に一様に分布する数を生成する.RandomVariateは指定された分布に対して定義された規則を使う.また,新しいメソッドと分布を定義するためのメカニズムも含まれている.

2引数インターフェースは,複数の乱数をすぐに得るのに便利である.さらに重要なことは,多数の擬似乱数を一度に非常に効率よく生成できるという利点がある.

0と1の間の10^7個の数の生成には数分の1秒しかかからない.
In[18]:=
Click for copyable input
Out[18]=
1回に1個ずつ10^7個の数を生成すると,ほぼ5倍の時間がかかる.
In[19]:=
Click for copyable input
Out[19]=

次元 から までの多次元配列では,要求される擬似乱数の合計 が生成されてから分割される.これによりできるだけ効率的に多次元配列が生成できる.それは,乱数値の合計数は可能な限り効率的に生成され,分割に必要な時間はほんのわずかしか必要ないからである.

100×100×100×10配列に必要な時間は,10^7個の数のベクトルに必要な時間にほぼ等しい.
In[20]:=
Click for copyable input
Out[20]=
In[21]:=
Click for copyable input
Out[21]=
同じ次元の配列を1度に10個の数で生成すると,数倍の時間がかかる.
In[22]:=
Click for copyable input
Out[22]=

統計分布では,一度に多くの数を生成することによる速度的な利点がさらに大きい.使用される同じ数の生成器から継承される効率面での利点に加え,統計分布の多くには,初等関数および特殊関数の評価のベクトル化からの利点もある.例えば,WeibullDistributionでは初等関数PowerTimesLogのベクトル評価が有利に働く.

10^5個のワイブル数の生成には,実質的に時間がかからない.
In[23]:=
Click for copyable input
Out[23]=
10^5個のワイブル数を一度に1つずつ生成すると,数秒かかる.
In[24]:=
Click for copyable input
Out[24]=

乱数生成は実地調査で便利なことがある.例えば,長い桁の数列の中で,ランダムな数列を探す場合等である.

これにより5つのランダムな10進数字を文字列に変換する.
In[25]:=
Click for copyable input
Out[25]=
次は の最初の100万桁を整数の文字列にする.
In[26]:=
Click for copyable input
これで, の最初の100万桁で先の10進数字が出現する場所が分かる.
In[27]:=
Click for copyable input
Out[27]=

乱数生成は,閉形式の結果が既知でない,あるいは計算的に難しいことが分かっている分布の推定におい手も非常に便利である.乱数行列の特性でその一例を示す.

下は,一様分布の実数の5×5行列が実数の固有値を持つ確率を推定する.
In[28]:=
Click for copyable input
Out[28]=
次は標準の正規分布の数の行列について同じことを行う.
In[29]:=
Click for copyable input
Out[29]=

多変量分布をシミュレートする例は,ベイズ統計で使用されるギブスサンプラーである[1].ギブスサンプラーは,他の座標を条件とする各座標の分布が既知であるならば,多変量分布からの値をシミュレートする方法を与える.ある制約条件下では,条件付き分布から繰返し抽出することにより構築された乱数ベクトルの分布が,本当の多変量分布に収束することがある.

次の例はCasellaとGeorgeの例[2]に対するギブスサンプラーを構築する.関心のある分布は2変量である. が与えられたときの の条件付き分布は二項分布, が与えられたときの の条件付き分布はベータ分布である.CasellaとGeorgeが述べているように,ギブスサンプラーを使った収束の検出および抽出のさまざまな方法が提案されている.簡単に言うと,収束は反復1000回以内に起ると仮定する.分布からの大きさ の標本は,1000回目の反復の後, 値とされる.しかし,この 値は従属値である.

二項およびベータの条件付き分布でサンプラーを定義する.
In[30]:=
Click for copyable input

ギブスサンプラーは,乱数生成の分布フレームワーク内部の分布オブジェクトとして定義することもできる.この特定の分布オブジェクトとしてのギブスサンプラーの例は,「分布生成器の定義」に記述してある.

は長さ10^4の標本である.
In[31]:=
Click for copyable input
次の棒グラフは,最初の次元の周辺分布を示す.
In[32]:=
Click for copyable input
Out[33]=
第2座標の周辺分布は,ヒストグラムで可視化できる.
In[34]:=
Click for copyable input
Out[34]=

条件付き分布は,それに対する十分なデータがある場合,仮定される二項およびベータ分布に近くマッチしなければならない.周辺分布の密度が高いときに最大量のデータが生じるので,それらの値は比較に使うことができる.次のグラフは,連続値の確率推定に幅.05のビンを使って,実際の条件付き分布と仮定された条件付き分布とを比較する.

の範囲での の実際の分布と理論的分布を比較する.
In[35]:=
Click for copyable input
Out[35]=
のときの の実際の分布と理論的分布を比較する.
In[36]:=
Click for copyable input
Out[36]=

任意精度の実数と複素数

RandomRealおよびRandomComplexは,デフォルトで機械精度数を生成する.RandomVariateはデフォルトでは連続分布に対する機械数を生成する.任意精度数はWorkingPrecisionオプションを設定して得ることができる.

オプション名
デフォルト値
WorkingPrecisionMachinePrecision計算で使用する算術の精度

RandomRealRandomComplexRandomVariateのオプション

このオプションは一様分布の実数および組込み分布からの複素数・実数に有効である.WorkingPrecisionはユーザ定義の分布に組み込むこともできる.

次は5から50までで精度25の実数である.
In[37]:=
Click for copyable input
Out[37]=
分布に従う精度50の数である.
In[38]:=
Click for copyable input
Out[38]=

WorkingPrecisionを増加すると,精度の損失が予想されたり非常に正確な結果が必要とされるシミュレーションに便利である.精度を増加すると,計算における精度の損失を推定することもできる.

区間での の計算における最悪の精度の損失を推定する.
In[39]:=
Click for copyable input
Out[39]=

入力の精度が指定されたWorkingPrecisionよりも小さい場合,関数により問題が警告される.入力の精度は,所望の精度の擬似乱数を生成するために人工的に増加される.

機械数7.5の精度は50より小さいので,警告が表示される.
In[40]:=
Click for copyable input
Out[40]=

WorkingPrecisionRandomIntegerのオプションではない.整数は無限精度を持つので,精度は関数名で完全に指定される.

WorkingPrecisionは整数の擬似乱数には意味をなさない.
In[41]:=
Click for copyable input
Out[41]=

ランダムな要素

RandomChoiceRandomSampleは,可能な要素のリストから擬似乱数的に選択する.要素は数値でも非数値でもよい.

RandomChoice[{e1,e2,}] のひとつを擬似乱数的に選択する
RandomChoice[list,n]list から n 個の要素を擬似乱数的に選択したリストを返す
RandomChoice[list,{n1,n2,}]list から擬似乱数的に選択した ××を返す
RandomChoice[{w1,w2,}->{e1,e2,}]
により重みの付いた擬似乱数的に選択した要素を返す
RandomChoice[wlist->elist,n]n 個の重み付きの擬似乱数的に選択した要素のリストを返す
RandomChoice[wlist->elist,{n1,n2,}]
重み付きの選択要素の ×× 配列を返す

リストからのランダムな選択

RandomSample[{e1,e2,},n] の中から n 個の無作為抽出を行う
RandomSample[{w1,w2,}->{e1,e2,},n]
重み を使って選ばれた の中から n 個の無作為抽出を行う
RandomSample[{e1,e2,}] の擬似ランダム置換を行う
RandomSample[wlist->elist]最初の重み wlist を使って elist の擬似ランダム置換を行う

リストからの無作為抽出

RandomChoiceRandomSampleの主な違いは,RandomChoice から復元抽出するのに対し,RandomSampleは非復元抽出するという点である.RandomChoiceで選ばれる要素の数は elist の要素の数で制限されず,要素 は複数回選択されてもよい.RandomSampleで返される標本の数は elist の要素の数で制限され,その標本の中で異なる要素が現れる回数は,elist でその要素が現れる回数で制限される.

RandomChoiceまたはRandomSampleの第1引数がリストならば,要素は同じ確率で選ばれる.重み指定は, の集合の分布を定義する.この重みは正でなければならないが,和が1である必要はない.重みでは,初期分布の の確率は である.RandomSampleは非復元抽出を行うので,各選択の後,残りの合計の重みに基づいて重みが内部的に更新される.

RandomChoiceは可能な結果の有限リストを持つ,独立で同一の分布の事象のシミュレーションに使うことができる.

公正なコイントス15回のシミュレーション.
In[42]:=
Click for copyable input
Out[42]=
5に対して負荷の付いた20回のさいころ投げ.
In[43]:=
Click for copyable input
Out[43]=

RandomChoiceは有限サポートを持つどの離散分布から観測を生成するためにも使うことができる.

以下はTriangularDistributionの離散分布から,ランダムな観測を生成する.
In[44]:=
Click for copyable input
Out[44]=
シミュレートされた1000個の点に対する経験的確率密度関数である.
In[45]:=
Click for copyable input
Out[46]=

RandomSampleは,結果のリストの各要素が1度だけしか観測されないような結果の有限集合からの観測をシミュレーションするために使うことができる.そのリストには異なる値が複数回出現してもよい.

次は青が80個,赤が45個入った容器から7回引くことをシミュレートする.
In[47]:=
Click for copyable input
Out[47]=

リストのすべての要素を無作為に抽出すると,ランダム置換になる.

次は1から10までの整数のランダム置換である.
In[48]:=
Click for copyable input
Out[48]=

要素に重みを割り当てると,置換において大きい重みを持つ値が,小さい値を持つ重みよりも先に現れる傾向にあるランダム置換となる.

データ値の平方で重みを付けたランダム置換.
In[49]:=
Click for copyable input
Out[49]=

重み付きあるいは重みなしの要素を持つ同じリストでは,RandomSample[#,1]&は分布がRandomChoiceに等しい.

大きさ1の10^5個の無作為抽出要素に対する経験的確率密度関数.
In[50]:=
Click for copyable input
In[51]:=
Click for copyable input
Out[51]=
In[52]:=
Click for copyable input
Out[52]=
分布的に等しいRandomChoiceに対する経験的分布である.
In[53]:=
Click for copyable input
Out[53]=

この2つの例に対する確率は互いに,また理論的値に非常に近い.

理論的確率である.
In[54]:=
Click for copyable input
Out[54]=

RandomSampleは,臨床試験等における,群への無作為割付けに使うこともできる.次では整数を使っているが,名前やID番号等の識別値を使うこともできる.

下は,20の元を同じ大きさの4つの群に無作為に割り付ける.
In[55]:=
Click for copyable input
Out[55]=

RandomChoiceおよびRandomSampleSeedRandomMethodオプションの変更の影響を受けることがある.組込みメソッドは「メソッド」で説明する.また,新しいメソッドを定義するメカニズムは「自分の生成器を定義する」で説明する.

種と局所化

擬似乱数生成では,明らかなランダム性のレベルを持つ数がアルゴリズムにより生成される.擬似乱数生成のためのメソッドでは,通常循環関係を使って現在の状態から数を生成し,次の数を生成する新しい状態を構築する.その状態は,アルゴリズムの循環関係を初期化するために使われる整数を生成器に撒くことで設定することができる.

「種」と呼ばれる初期値が与えられると,擬似乱数生成器は完全に確定的になる.多くの場合,乱数生成に対する種を局所的あるいは大域的に設定し,「乱数」値の数列を継続的に得ることが望ましい.大域的に設定すると,新しい種を明示的に設定しない限り,この種が今後の擬似乱数に影響する.局所的に設定すると,種は局所化されたコード内の乱数と要素の生成にのみ影響を与える.

BlockRandom[expr]すべての擬似乱数生成を局所化して expr を評価する
SeedRandom[n]n を種として使い,擬似乱数生成をリセットする
SeedRandom[]時間と現在のWolfram言語セッションのある属性を種として使い,生成器をリセットする

局所化と種関数

SeedRandom関数は乱数生成器に種を蒔く方法を提供する.SeedRandomだけを使うと,乱数生成器用の種は大域的に設定される.BlockRandom関数は,大域的な状態に影響を与えずに乱数生成器用の種を局所的に設定あるいは変更する方法を提供する.

次は大域的に乱数生成器の種を蒔く.
In[56]:=
Click for copyable input
Out[56]=

次では最初のRandomRealBlockRandom内で生成され,2つ目がBlockRandomの外側で生成されているため,2つの異なる数を返す.

2つ目のRandomRealは種を使って生成されない.
In[57]:=
Click for copyable input
Out[57]=

SeedRandomは乱数生成を切り替えるメカニズムも提供する.

オプション名
デフォルト値
MethodAutomatic蒔いて使用するメソッド

SeedRandomのオプション

個々の生成器にはMethodオプションでその生成器を指定することにより,直接種を蒔くことができる.すべての生成器には,設定Method->Allで種を蒔くことができる.

ここではデフォルトの生成器の種はであるが,生成器には種がない.
In[58]:=
Click for copyable input
Out[58]=
に種を蒔くと,異なる乱数が生成される.
In[59]:=
Click for copyable input
Out[59]=

並列計算のSeedRandomおよびBlockRandom

並列計算でSeedRandomおよびBlockRandomというコマンドを使うと少し異なる点がある.並列計算内部では,これらのコマンドは現行のスレッドで使用されている生成器にのみ影響を及ぼす.通常,これらのコマンドは並列計算の前,あるいは並列計算を包み込んで使った方がよい.

並列計算では,各スレッドで乱数を生成する生成器は他のスレッドの生成器に非依存である方が有利である.Wolframシステムでは,並列計算で使われる各スレッドに対して,0から始まり通常は順に$ProcessorCountまで至る一意の指標が与えられている.これは各スレッドに異なるシードや生成器を与えるために使われる.

次の表はこれらが逐次計算と並列計算で使われる場合の差分を説明したものである.

コマンド
逐次計算
並列計算
SeedRandom[seed]現行の逐次型乱数生成器に seed という種を並列型生成器に seed + i という種を播く.ここで i は並列スレッドの指標である現行スレッドの乱数生成器にのみ seed という種を播く
SeedRandom[seed,Method->"ParallelGenerator"]並列型生成期に seed + i という種を播く.ここで i は並列スレッドの指標である影響なし
SeedRandom[Method->method]逐次型の乱数生成器に対するメソッドから method へと変更する現行のスレッドに対する乱数生成器のみに対するメソッドから method へと変更する
BlockRandom[expr]局所化されたすべての擬似乱数生成器で expr を評価する局所化された現行のスレッドに対する擬似乱数生成器のみで expr を評価する

逐次計算および並列計算におけるSeedRandomおよびBlockRandom

以下でCompiledFunctionを定義する.これは,リストが与えられたときに並列で実行する 個のサンプルを使って,1/4円の面積を近似する.
In[60]:=
Click for copyable input
次は,すべての生成器の種を播いた後,CompiledFunctionを並列で実行する.
In[61]:=
Click for copyable input
Out[61]=
再び実行するが,今回は並列計算がBlockRandom内部で行われる.
In[62]:=
Click for copyable input
Out[62]=

同じ種であるにもかかわらず結果が異なっている.この差分のほとんどは順序付けにある.計算が繰り返される場合,並列スケジューラは1つのスレッドを実行する前に別のスレッドを実行することがあるからである.

結果を比較する.
In[63]:=
Click for copyable input
Out[63]=

同じ結果のすべてまでではないが,多くが両方の結果に現れている.これは,計算が繰り返されるとき,指定されたスレッドが正確に同じ回数だけ使われるとは限らないためである.

前の並列計算はBlockRandomの内部で行われたが,並列ジェネレータは以前の状態に戻されているため,再び実行することは実質的に繰り返すことと同じである.
In[64]:=
Click for copyable input
Out[65]=

並列計算の内部でSeedRandomおよびBlockRandomを使う場合は注意しなければならない.同じスレッドで行われた異なる部分が,結果として同一の結果になる可能性があるからである.

以下では,先ほどと同じCompiledFunctionを定義しているが,RandomRealを使う前にSeedRandomを実行するようになっている.
In[66]:=
Click for copyable input
並列でCompiledFunctionを実行する.
In[67]:=
Click for copyable input
Out[80]=

結果の中のいくつかが同じになっているのが分かる.これはUnionを使って調べることができる.

以下で,結果から一つずつしかない和を得る.
In[81]:=
Click for copyable input
Out[81]=

この場合,20個のうちで単独の和は8個しかない.もしこれを実行したら,和集合の長さは通常,使用中のマシンに搭載のプロセッサの数に等しくなる.各プロセッサの生成器には,使用される前に毎回種が播かれるため,またそれぞれの場合のRandomRealの使用方法が同じであるため,結果が同一になるのである.

前と同じようにCompiledFunctionを定義するが,今回はBlockRandomの内部でRandomRealを使う.
In[91]:=
Click for copyable input
並列でCompiledFunctionを実行する.
In[92]:=
Click for copyable input
Out[92]=
結果から一つずつしかない和を得る.
In[93]:=
Click for copyable input
Out[93]=

並列計算の内部でSeedRandomを使ってできることのひとつに,生成器の設定がある.各スレッドの生成器が,デフォルトの生成器で異なる種を持つように指定したいとする.

次で乱数生成器をメソッドに変更し,それに種 s を播くコンパイルされた関数を定義する.
In[6]:=
Click for copyable input
それぞれの生成器に対してランダムに選ばれた種を与える.
In[2]:=
Click for copyable input
並列でCompiledFunctionを実行する.並列の乱数生成器だけがこの影響を受ける.
In[3]:=
Click for copyable input
並列で面積近似関数を実行するときは次の生成器が使われる.
In[9]:=
Click for copyable input
Out[9]=

You can verify that these generators were used by comparing to a serial computation where the generator is set the same way.

並列計算で設定されたのと同じように局所的に生成器を設定して,逐次計算する.
In[10]:=
Click for copyable input
Out[10]=

The parallel result is just a permutation of this.

並列計算の結果が逐次計算の結果の置換であることを検証する.
In[29]:=
Click for copyable input
Out[31]=
In[32]:=
Click for copyable input
Out[32]=

同じ生成器で種を替えるだけでは生成された数に何らかの相関関係が発生するとは限らないため,このように生成器を設定することはお勧めできない.

並列の生成器を設定するためのより簡単で信頼できる方法は,「メソッド」セクションで説明する"ParallelGenerator"メソッドで提供されている.

メソッド

すべてのシステムで,5つの擬似乱数生成メソッドが使用できる.これら5つのうち,メルセンヌ・ツイスタ(Mersenne Twister)法は逐次・並列の両方で提供されている.6番目のプラットフォーム依存のメソッドはIntelベースのシステムで利用できる.メソッドの名前は並列計算の生成器の取扱いのために使われる.「自分の生成器を定義する」で説明してある新しいメソッド定義のためのフレームワークも含まれている.

"Congruential"線形合同法(低質のランダム性)
"ExtendedCA"拡張セルオートマトン法(デフォルト)
"Legacy"Wolframシステム 6.0より前のデフォルト法
"MersenneTwister"メルセンヌ・ツイスタシフトレジスタ法
"MKL"インテルMKL(インテルベースのシステム)
"ParallelGenerator"並列計算のために生成器を初期化し種を播くために使われる
"ParallelMersenneTwister"周期の1024個のメルセンヌ・ツイスタ生成器の集合
"Rule30CA"ウルフラム則30生成法

組込みメソッド

シードが2020のときの各メソッドからの整数の擬似乱数.
In[54]:=
Click for copyable input
Out[54]=
同じシードから実数の擬似乱数を与える.
In[55]:=
Click for copyable input
Out[55]=

Congruential

は線形合同法を使う.これは擬似乱数が から得られる0と1である,最も単純な擬似乱数生成法のひとつである. は下の漸化式

であり,固定された整数 はそれぞれ乗数,増分,法と呼ばれる.増分が0ならば,生成器は乗法合同法である. の値はオプションを使ってメソッドに設定することができる.

オプション名
デフォルト値
"Bits"Automaticビットから構成される数に使用するビットの範囲を指定する
"Multiplier"1283839219676404755乗数の値
"Increment"0増分の値
"Modulus"2305843009213693951法の値
"ConvertToRealsDirectly"True合同関係から直接実数を構築するかどうか

Method のオプション

線形合同法は周期的であり,多数の乱数値が必要な場合は特にランダム性の質が低下する傾向にある.実数が合同関係から直接生成される場合,周期は 以下である.

デフォルトのオプション値は,大きい周期を持ち,64ビット効率に対応するよう選ばれる.デフォルトオプションでは,生成器は,合同数生成に関する継承問題があるにもかかわらず,乱数性をみる多数の標準テストに合格する.

乗法合同法から40個の数を生成する.
In[56]:=
Click for copyable input

乗法合同法の周期は,法に互いに素である法と等しいかそれ以下の正の整数の数により制限される.この上限は,オイラーのトーティエント関数である.

法が63では,周期は最大で36である.
In[57]:=
Click for copyable input
Out[57]=

実際の周期は を法とする となるような最小の整数 を見付けることで決めることができる.

乗数11,法63のときの周期は6である.
In[58]:=
Click for copyable input
Out[58]=
データを6つの要素の集合に分割すると再帰が分かる.
In[59]:=
Click for copyable input
Out[59]=

生成された数列をプロットすると,異なる数字がグラフィカルに表示される.

合同法からの1000個の値をプロットする.
In[60]:=
Click for copyable input
Out[60]=

Falseに設定されると,数列の要素から一度に8ビット取り52ビットの機械精度数を構築することにより実数が生成される.この方法で生成された合同数はまだ循環するが,循環は初期の合同関係の反復ではなくビットパターンの反復に依存する.

オプションはAutomatic,非零の整数,ビットからの数を構築するために使われる法 のビット範囲を指定する2つの非零の整数のリストのいずれかを取る.Automaticではが使われるが, が2のベキ乗のときはが使われる.

ExtendedCA

デフォルトのメソッドは,セルオートマトンを利用して高質の擬似乱数を生成する.この方法は,特殊な5近傍則を使うので,それぞれの新規セルは,直前のステップの5つの非隣接セルに依存する.

セルオートマトンベースの乱数生成は,決定的規則に従って0と1の状態ベクトルを進化させる.指定されたセルオートマトンでは,新しい状態ベクトルの指定の位置にある要素(またはセル)は,古い状態ベクトルにおけるそのセルの一定の隣接セルにより決められる.状態ベクトルのセルの部分集合は,ランダムなビットとして出力され,それから擬似乱数が生成される.

で使用されるセルオートマトンは,非常に高レベルのランダム性を生成する.そのランダム性は非常に高度なので,出力ですべての単独セルを使った場合でさえも,1つのセルと5つの前のセルの間に明らかな相関があるにもかかわらず,多数のランダム性テストにパスするようなビットのストリームを生成する.

状態ベクトルの大きさ,飛ばされるセル,開始セルを変更するためのオプションが含まれている.デフォルトは質と速さで選ばれ,通常このオプションを変更する必要はない.

オプション名
デフォルト値
"Size"8064の乗数としての状態ベクトルの大きさ
"Skip"4飛ばすセルの数
"Start"0どのセルから始めるか

Method のオプション

使用される状態ベクトルの長さはデフォルトでセルに設定されている.64の乗数がオプションで指定できる.5近傍ルールを使ってセルオートマトンを進化させることにより状態ベクトルが計算されると,ビットから乱数に対するビットが選ばれる.

実際には,各状態ベクトルで4つおきにセルを使うことは,非常に厳しいランダム性テストにパスするのに十分であると分かっている.これはオプションに使われるデフォルトである.さらに高速な乱数生成法では,設定が2もしくは1にもなるが,乱数の質は低下する.

より大きいおよびと互角であるオプションは,並列計算で使える非依存の生成器族を設定するのに便利である.

はデフォルトの数生成法である.
In[61]:=
Click for copyable input
Out[61]=
In[62]:=
Click for copyable input
Out[62]=

Legacy

メソッドは,Mathematica バージョン6.0より前ではRandomにより呼び出される生成法を使用する.実数にはマーサグリア・ザマンボロー付き減算ジェネレータが使われる.整数の生成法にはウルフラム則30のセルオートマトンジェネレータを使う.規則30ジェネレータは,小さい整数の場合は直接,大きい整数の場合は一定のビットを生成するために使われる.

これはメソッドで得られたRandomRealRandomIntegerの値である.
In[63]:=
Click for copyable input
Out[63]=
同等のRandom呼出しでも,同じ値が得られる.
In[64]:=
Click for copyable input
Out[64]=

バージョン6.0以前に生成された列との一貫性を保障するため,Automaticメソッドの種の集合がメソッドにも適用される.

メソッドにはオプションがない.

MersenneTwister

はMatsumotoとNishimura [3][4]によるメルセンヌツイスタ生成器を使用する.メルセンヌツイスタは,周期の一般化されたフィードバックシフトレジスタ生成器である.

メルセンヌツイスタ生成器から5つの乱数を与える.
In[65]:=
Click for copyable input
Out[65]=

メソッドにはオプションがない.

MKL

メソッドは,IntelのMKLライブラリで提供されている乱数生成器を使う.MKLライブラリは,プラットフォーム依存である.メソッドは,Microsoft Windows(32ビット,64ビット),Linux x86(32ビット,64ビット),Linux Itaniumの各システムで使用可能である.

オプション名
デフォルト値
MethodAutomatic使用するMKL生成器

Method のオプション

"MCG31"31ビット乗算合同法
"MCG59"59ビット乗算合同法
"MRG32K3A"三次の2つの要素を持つ結合多重再帰発生器
"MersenneTwister"メルセンヌツイスタシフトレジスタ生成器
"R250"一般化フィードバックシフトレジスタ生成器
"WichmannHill"WichmannHillの結合乗算合同法
"Niederreiter"Niederreiterの低食い違い量列
"Sobol"Sobolの低食い違い量列

メソッド

最初の6つの方法は一様乱数生成法である.はNiederreiter列とSobol列を生成する.これらの数列は非一様であり,数値法で便利な場合もある内在構造を持つ.例えば,これらの数列は通常多次元のモンテカルロ積分で速く収束する.

以下は二次元のNiederreiter列を示す.
In[66]:=
Click for copyable input
Out[66]=
二次元のSobol列の構造を示す.
In[67]:=
Click for copyable input
Out[67]=

Rule30CA

メソッドはウルフラム則30のセルオートマトン生成器を使う.以下の関係

を使い,0と1の状態ベクトルを進化させることによりビットを得る.ここで は時間 におけるセル の値である.

オプション名
デフォルト値
"Size"929の乗数としての状態ベクトルの大きさ

Method のオプション

使用される状態ベクトルの長さは,デフォルトでセルに設定されている.29の乗数はオプションで指定することができる.

を使って整数の乱数の2x3x4テンソルを与える.
In[68]:=
Click for copyable input
Out[68]=

メソッドは,各状態ベクトルから最初のビットしか使わないので,各状態ベクトルから複数のビットを使用するメソッドよりも遅い.

ParallelMersenneTwister

ではMatsumotoとNishimura [3][4]による,「Dynamic Creator」プログラムdcmt[19]を使って選ばれたパラメータを持つメルセンヌ・ツイスタ生成器が一組使われる.このプログラムは,互いに素であるために非依存の結果を出すメルセンヌ・ツイスタ生成器に対するパラメータを計算する.そのパラメータは,周期でメルセンヌ・ツイスタの一般フィードバックシフトレジスタ生成器を作成するために計算される.

どの生成器を使うかを決めるオプションが含まれている.

オプション名
デフォルト値
"Index"00から1023間でのどの生成器を使うか

Method のオプション

次は,異なる並列メルセンヌ・ツイスタ生成器から2500個の乱数を含む2つの集合を与え,そのペアのプロットを点で作成する.
In[9]:=
Click for copyable input
Out[10]=

この2つの生成器で生成された乱数には明確な相関関係はない.相関関係がない上,スピードも遅いため,この生成期の集合が並列計算のためのデフォルトの生成器となっている.

ParallelGenerator

は並列計算に使われる生成器のスピード調整や変更を可能にするコントローラメソッドである.

どの生成器集合を使うかを決めるためのオプションが含まれている.

オプション名
デフォルト値
MethodAutomaticどの非依存生成器を使うか

Method のオプション

メソッドに与えることのできるMethodオプションの値は,組込みのパラメータ化されたオプションか非負の整数に対してランダムな生成器指定を与える関数を指定する文字列である.並列計算で使用されるそれぞれのスレッドには,ゼロから始まり通常$ProcessorCountまで続く一意の指標が与えられる.この指標はそれぞれのスレッドに異なる種と生成器を与えるために使われる.

"ParallelMersenneTwister"周期の並列メルセンヌ・ツイスタ法
"ExtendedCA"異なる開始位置の拡張セルオートマトン法
fi番目のスレッドに使われる生成器
"Default"デフォルトのメソッドに戻す

並列生成器のメソッド

独立した2セットの高質生成器を得る便利な方法として,文字列のショートカットが提供されている.

を使うことは,関数f=Function[{i},{"ParallelMersenneTwister","Index"->i}]を使うことに等しい.生成器が高速であり,高品質の乱数を生成するため,これが並列計算のデフォルトとなっている.

を使うことは,以下で定義した関数 f を使用中のマシンのプロセッサの数で使うことと通常等しい.

はメソッドをデフォルトのメソッドに再設定する.

次でのオプションMethod->"ExtendedCA"でのデフォルト関数を定義する.
In[11]:=
Click for copyable input
Out[11]=

パラメータは,マシン上のすべての$ProcessorCountのプロセッサを使っても,デフォルトの逐次乱数生成器と同じくらいよい乱数が得られるように選ばれる.

メソッドも少し異なる方法で生成器の種蒔きを行う.各プロセッサに同じ種を使う代りに,SeedRandom[seed, Method->"ParallelGenerator"]は各スレッドに対して i はそのスレッドについての一意の指標)を使う.これにより,たとえ各スレッドの生成器を同じに設定しても(例:Method->Function[{i},"ExtendedCA"]),異なるスレッドから異なる乱数が得られる.しかし,異なる種であるとしても乱数が思いもよらない相関性を持つことがあるため,この方法はお勧めできない.

一般に,異なるスレッドに対する生成器メソッドを与える関数 f は,正当な乱数生成法を返す.

リストを与えると並列で実行するCompiledFunctionである.
In[1]:=
Click for copyable input
並列生成器の種を播く.
In[2]:=
Click for copyable input
CompiledFunctionを並列で実行する.
In[3]:=
Click for copyable input
Out[3]=
次は0から7の間の指標に対して異なる生成法を与える関数を定義する.
In[4]:=
Click for copyable input
次は並列生成器を,関数によって与えられた生成器に変更し,それに種を播く.
In[12]:=
Click for copyable input
選ばれた生成器を使って,コンパイルされた関数を並列で実行する.
In[13]:=
Click for copyable input
Out[13]=
次は逐次に計算を実行し,局所的に生成器を関数によって与えられたものに設定する.
In[14]:=
Click for copyable input
Out[14]=
結果は.順序を別にすれば同じである.
In[15]:=
Click for copyable input
Out[15]=

並列生成器をデフォルトメソッドに戻すためには,明示的にメソッドオプションを与える必要がある.そうしないと,種を変更するだけになってしまう.

次で並列生成器をデフォルトのメソッドに戻す.
In[16]:=
Click for copyable input

自分の生成器を定義する

メソッドは正しいテンプレートに従っている限り,乱数フレームワークにプラグインすることができる.生成期のオブジェクトは という形式である.ここで gsym は生成器を同定し,規則が接続されているシンボルである.data は実質的に,生成器の定義に関連付けられたトップレベルの評価に対してプライベートである.

生成器の初期化はの呼出しにより実行される.

Random`InitializeGenerator[gsym,opts]
生成器 gsym をオプション opts で初期化する

生成器初期化関数

は形式 の生成器オブジェクト gobj を返すことが想定されている.

生成器はランダムなビットストリーム,ランダムな整数,ランダムな実数の生成をサポートすることができる.生成器がビットストリームをサポートする場合,実数と整数はビットストリームの変換により生成できる.メソッドの設定時には,何をどのようにサポートするかを決定するための属性が問い合される.

GeneratesBitsQメソッドがビットを生成する場合はTrueに設定
GeneratesIntegersQメソッドが指定された範囲の整数を生成するならばTrueに設定
GeneratesRealsQメソッドが指定された範囲と精度で実数を生成するならばTrueに設定

生成器の属性

ビットストリームがサポートされる場合,n 個のランダムなビット,あるいは項目が0か1の長さ nbits のリストを返すことが想定される.

ランダムな整数がサポートされる場合,は範囲n 個のランダムな整数のリストを返すことが想定される.結果が範囲を超えると警告メッセージが出される.

ランダムな実数がサポートされる場合,は範囲で精度 precn 個の実数のリストを返すことが想定される.結果が範囲を超えたり,異なる精度であったりした場合は,警告メッセージが出される.

生成関数のすべてに対して,返されるものはである.res は正しいタイプの結果であり,gobj は新しい生成器オブジェクト(状態の変化を反映して)である.

整数 seed に対する種まきは により実行される.は新しい生成器オブジェクトを返すものと想定される.

例題:乗算合同法

次の例では,乗算合同法を定義する.乗算合同法は,漸化式

に従う.以下で定義されるように,生成器は実数の生成だけができる.

生成器のデフォルトオプションを設定する.
In[69]:=
Click for copyable input

生成器の初期化は乗数と法の値を抽出する.この値のいずれかが正の整数でない場合,初期化は失敗する.

次で生成器を初期化する.
In[70]:=
Click for copyable input

カーネルからの呼出しは,実質的にCatchにラップされている.問題が生じた場合は,初期化コードにThrowを使って簡単に終了することができる.

が実数を生成するよう設定する.
In[71]:=
Click for copyable input
漸化式を使い,生成器に種を蒔く.
In[72]:=
Click for copyable input

実数生成器は所望の数の実数と新しい生成器を返す.新しい生成器の種は,漸化式に基づき更新される.

これで実数生成器を定義する.
In[73]:=
Click for copyable input
生成器を使い,10個の実数を生成する.
In[74]:=
Click for copyable input
Out[74]=
この生成器は整数に対しては定義されていない.
In[75]:=
Click for copyable input
Out[75]=

例題:BlumBlumShub法

Blum Blum Shub法は暗号化目的の擬似乱数ビットを生成するための二次合同法である[5].合同式は指定された素数 に対する法 × である.

生成器のデフォルトオプションを指定する.
In[76]:=
Click for copyable input
以下で補助関数と生成器に対するエラーメッセージを定義する.
In[77]:=
Click for copyable input
In[78]:=
Click for copyable input
In[79]:=
Click for copyable input
In[80]:=
Click for copyable input

生成器の初期化でオプション値が抽出され,実際の生成器を呼び出す前に,必要に応じてエラーメッセージが出される.

次で生成器を初期化する.
In[81]:=
Click for copyable input
がビット生成器であることを指定し,ビット幅を決める.
In[82]:=
Click for copyable input
In[83]:=
Click for copyable input
次で生成器に種を蒔く.
In[84]:=
Click for copyable input
ビット生成器を定義する.
In[85]:=
Click for copyable input
生成器を使い,5つの整数と5つの実数を生成する.
In[86]:=
Click for copyable input
Out[86]=

統計分布

非一様分布からランダムな変量を生成する目的は,0と1の間のランダムな一様変量を生成し,所望の分布の乱数値の累積分布関数の逆関数を計算することにある.しかし,実際には,直接これに従うと多数の乱数変量が必要な場合に計算的に非常に厳しくなる.特に累積分布関数の逆関数が複雑で,閉形式で表せない場合はそうである.

そのような場合,累積分布関数の直接の逆関数にはテーブル検索,分布関係に基づく直接構造,棄却採択法の方が効率的であることが多い.あるレベルにおいては,このような方法はすべて一様分布に従うRandomReal値,非一様分布のRandomInteger値,重み付きのRandomChoiceの観測値,あるいはこれらすべての値の組合せに依存する.その結果,SeedRandomで設定された方法は,統計分布の乱数観測値に影響を及ぼす.

すべての組込みの統計分布からの乱数観測値は,RandomVariateを使って生成することができる.Wolfram言語の多くの分布に対してRandomVariateで使われるメソッドは,Gentle [6]あるいは他の文献で提唱もしくは記述されているメソッドに従う.

RandomVariate[dist]連続分布 dist からの乱数を与える
RandomVariate[dist,n]dist からの実数の擬似乱数 n 個のリストを与える
RandomVariate[dist,{n1,n2,}]dist からの実数の擬似乱数の ××配列を与える

統計分布からの乱数値の生成

統計分布からの観測値はRandomVariateで得ることができる.これには一変量および多変量分布,連続および離散分布,パラメトリックおよび導出分布,データにより定義された分布を含むすべての組込み分布や構文が含まれる.

次で連続分布および離散分布の数値を生成する.
In[2]:=
Click for copyable input
Out[2]=

WorkingPrecisionは,範囲上の一様乱数に対するのと同様に,連続分布に対するより精度の高い値を得るために使われる.

精度30のベータ分布変量である.
In[3]:=
Click for copyable input
Out[3]=

多変量分布からの乱数値も同様に生成することができる.

2変量正規分布からの乱数ベクトルである.
In[90]:=
Click for copyable input
Out[90]=
多変量正規分布からの乱数ベクトルである.
In[91]:=
Click for copyable input
Out[91]=
ここでは,乱数値は確率密度分布により定義された分布から生成される.
In[8]:=
Click for copyable input
Out[9]=

このセクションでは,確率変数生成のためのメソッドについて,そのようなメソッドがWolfram言語のどこで使われるかに関していくつかの具体的な例を交えて説明する.

連続分布

累積密度関数の逆関数に初等関数しか含まれていないような一変量分布では,一様乱数に対する累積密度関数の逆関数の直接計算が通常使われる.これは一様分布に従う確率変数からの直接構築と見ることができる.これに分類される連続分布の中にはCauchyDistributionExponentialDistributionExtremeValueDistributionGumbelDistributionLaplaceDistributionLogisticDistributionParetoDistributionRayleighDistributionTriangularDistributionWeibullDistributionがある.

複数の一様変数から,あるいは一様分布以外の変数からの単独の確率変数の直接構築も使用される.正規変数は,BoxMüller法を使い一様乱数のペアから,ペアとして生成される.例えば,HalfNormalDistributionLogNormalDistributionMultinormalDistributionの確率変数は,正規変数の直接変換により得られる.

InverseGaussianDistributionは,正規変数と一様変数を含むacceptancecomplement法を使用する.この方法はMichael,Schucany,Haasによるものであり,Gentle [6]に記述されている. MaxwellDistribution変数はChiDistribution変数から構築される.カイ変数自体は,GammaDistributionの特殊形であるChiSquareDistributionから得られる.

ほとんどの場合,FRatioDistributionは単独のベータ変数乱数からそれぞれの乱数値を構築する.自由度が小さい場合,FRatioDistribution変数は,ベータ構築で発生する可能性のある0による割り算の可能性を避けるために,ガンマ変数のペアから生成される.

NoncentralChiSquareDistribution[ν,λ],の変数生成では,分布の加法的性質を使い,非整数 に対する累積密度関数の高価な逆関数の計算を避ける.この加法的性質は,例えば Johnson,Kotz,Blakrishnan [7]で与えられる.のとき,非心 変数は平均,分散1の正規変数の平方として生成することができる. では,非心 変数は中心および非心 確率変数の和として得られる.のとき, はX~およびY~ならば 分布に従う.この関係は のときは使えない.その場合,構文は および である となる.ここで が0に近付くときの極限分布の非心 分布である.極限分布 はポワソン変数と 変数を合わせたものであり,0において非零の確率質量を,正の値の時は連続密度を持つ.NoncentralFRatioDistribution変数は,中心および非心の 変数1つずつから得られる.

多変量統計パッケージのWishartDistributionでは,行列はSmithとHockingの方法[8]で生成される.この方法は,カイ分布に従った対角要素と正規分布に従った非対角要素を持つ行列からウィッシャート行列を構築する.

NoncentralStudentTDistributionHotellingTSquareDistributionMultivariateTDistributionはそれぞれ一変量確率変数からの直接構築を使う.

GammaDistributionBetaDistributionStudentTDistributionはある程度まで棄却採択法を使う.

GammaDistribution[α,β]指数変数は のときに生成される.それ以外のときは,Cheng,Feast [9]によるメソッドおよびAhrens,Dieter [10]によるメソッドが使われる.

ベータ変数はベータパラメータ および の値により複数の方法を切り替えることにより構築される.両方のパラメータが1ならば,一様確率変数が生成される.ベータパラメータのうちのひとつが1ならば,累積密度関数の閉形式の逆関数評価が使われる.それ以外のとき,RandomVariateはJöhnk [11],Cheng [12],Atkinson [13]による棄却採択法のいずれかを使う.2つのガンマからの構築よりも棄却採択法を使うことの利点の一例を下に挙げる.直接棄却採択法はガンマペア構築法のほぼ2倍の速さである.

ベータ変数の直接構築と棄却採択法とを比べる.
In[92]:=
Click for copyable input
Out[92]=
In[93]:=
Click for copyable input
Out[93]=

StudentTDistributionでは,RandomVariateで使用されるメソッドはBailey [14]のPolar Rejection法である.この方法は以下で示す通り,正規変数および 変数からの直接構築よりも効率的である.直接構築は100万のスチューデントの 変数に対してPolar法のほぼ1.5倍の時間がかかる.

スチューデントの に対する直接構築とBaileyのPolar Rejection法を比較する.
In[94]:=
Click for copyable input
Out[94]=
In[95]:=
Click for copyable input
Out[95]=

離散分布

GeometricDistributionBetaBinomialDistributionBetaNegativeBinomialDistributionはすべて直接構築を使う.GeometricDistribution変数はとして生成され,UniformDistribution[0,1]に従う.BetaBinomialDistributionおよびBetaNegativeBinomialDistributionは,確率パラメータをBetaDistributionの確率変数とするBinomialDistribution変数およびNegativeBinomialDistribution変数から構築される.

整数の乱数生成で表検索が使われるときは,重みに対する分布の確率質量関数を使い,RandomChoiceにより実装される.ほとんどの離散分布は,所望の標本サイズが与えられたときに,重みのリストの構築が効果になると見込まれると,別の方法に切り替える.例えば, が1に近付くと,LogSeriesDistribution[p]は直接構築に切り替える.ここで, および は区間 [15]上で一様分布に従う.また,パラメータと標本の大きさによっては,NegativeBinomialDistribution[n,p]はポワソンとガンマを合わせた構築法に切り替える.これは平均がガンマ分布に従うポワソン変数[6]である.

BinomialDistributionHypergeometricDistributionPoissonDistributionはすべて,確率分布関数値の計算のオーバーヘッドが所望の乱数値の数と比較して小さければ,密度関数からの直接抽出に依存する.確率分布関数値の直接計算においてオーバーフローやアンダーフローが起きると,棄却採択法も変数の生成が行えるので,乱数が生成できるパラメータ値の範囲が拡張される.

二項分布および超幾何分布は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[μ]NormalDistribution[μ,]に傾くことを利用している.

ZipfDistributionの乱数値は,Devroye [15]により記述され手いる棄却採択法で生成される.このメソッドは一様変数のペアおよび基礎演算以外ではFloorと非整数ベキしか含まない検定を使い,ジップ分布値を効率的に得る.

分布生成器の定義

Wolfram言語には多くの分布構文が含まれているため,他の分布と同様に扱うことのできる新しい分布オブジェクトを定義することが可能である.これには乱数生成も含まれる.しかし,ある分布から乱数を生成するだけでよく,そのためのアルゴリズムもあるとしたら,乱数生成器だけを定義した方がよい場合もある.そのような分布生成器の定義は,の規則でサポートされている.

Random`DistributionVector[expr,n,prec]
精度 precdist から n 個の観測値を生成するための規則を定義する

分布から乱数生成を定義するための関数

は指定された精度で指定された長さのベクトルを返す.式 expr は完全に定義された分布ではないので,数値はRandomVariateではなくRandomRealまたはRandomIntegerで生成される.精度がInfinityの場合は,数値はRandomIntegerで生成される.それ以外の場合はRandomRealで生成される.

関数から乱数値を生成するための規則は,一般に分布の頭部のTagSetで定義される.分布自体にパラメータが含まれることもある.簡単な例として,以下で区間上の一様分布を表すに対する規則を定義する.

In[96]:=
Click for copyable input

これでの乱数が,RandomRealで生成できるようになった.

以下はから機械精度数および20桁精度の数を与える.
In[97]:=
Click for copyable input
Out[97]=

行列および高次元のテンソルもRandomRealで直接生成することができる.RandomRealに与えられる定義を使い,所望の乱数値の合計数を生成し,それを指定された次元に分割する.

数の3x4行列である.
In[98]:=
Click for copyable input
Out[98]=

離散分布生成器も同様に定義することができる.大きく異なる点は,の精度引数がInfinityであるということである.の離散型で,簡単な例を示す.

In[99]:=
Click for copyable input

の乱数値は,RandomIntegerから得られるようになった.

の10個の数である.
In[100]:=
Click for copyable input
Out[100]=

前の例では分布生成器の定義に対する基本的なフレームワークを示したが,分布自体は特に面白いものではない.実際,この2つの場合では新しい分布のシンボルに定義を加える代りに,RandomVariateから値だけを生成し,最終結果にを掛けた方が簡単である.次の例ではシンボルに定義を加えた方が便利であるような,少し複雑な分布を示す.

例題:逆関数による正規分布

一般の一変量分布から乱数値を生成する教科書的な定義には下の2つのステップが含まれる.

  • 区間における一様乱数 を生成する
    • 関心のある分布に対する の累積分布関数の逆関数を計算する

    このプロセスを例示するために,ここでは正規分布を使う.この方法を使ってランダムな正規変数を生成するためには,点 において正規分布に対するQuantileから開始し, を0と1の間の一様乱数で置き換えるとよい.

    平均 ,標準偏差 の正規分布に対するQuantile関数.
    In[101]:=
    Click for copyable input
    Out[101]=

    これで新しい分布オブジェクトを使い,累積分布関数の逆関数を使う正規乱数生成器を定義することができる.

    逆関数により正規分布の生成を定義する.
    In[102]:=
    Click for copyable input
    次は逆関数により生成された10個の正規乱数である.
    In[103]:=
    Click for copyable input
    Out[103]=
    標本平均および標準偏差を伴う10^4個の正規乱数の標本である.
    In[104]:=
    Click for copyable input
    Out[104]=
    In[105]:=
    Click for copyable input
    Out[105]=

    正規分布は,直接倒置以外の方法が好まれる場合のひとつである.累積分布関数の逆関数が擬似正規乱数生成に対して完全に有効なメソッドである反面,完全に効率的であるわけではない.InverseErfの数値評価は,NormalDistributionによって使用されるBoxMüller法で要求される正弦関数や対数関数の評価よりも計算的に効果である.

    組込みの方法では同じ数の値に対してほとんど時間がかからない.
    In[106]:=
    Click for copyable input
    Out[106]=
    In[107]:=
    Click for copyable input
    Out[107]=
    In[108]:=
    Click for copyable input

例題:円板上の一様分布

は多次元分布の生成器を定義するために使うこともできる.例えば,単位円板上の一様分布からの乱数点があるとすると,である実数点の集合が望まれる.そのような乱数点は次のように構築することができる.

  • で一様分布するランダムな角度 を生成する
  • で一様分布するランダムベクトル を生成する
    • を返す

    半径 の円板上で一様分布する点を生成するために,戻った順序対に を掛けることができる.

    次は半径 の一様円板に対する生成器を定義する.
    In[109]:=
    Click for copyable input
    半径2の円板からのランダムな順序対である.
    In[110]:=
    Click for copyable input
    Out[110]=
    次は,この円板上の10^4個の生成点の分布を可視化する.
    In[111]:=
    Click for copyable input
    Out[111]=

例題:ギブスサンプラー

ギブスサンプラーは分布生成器として定義することもできる.例として,ベータ分布と二項分布を組み合せたギブスサンプラーを見てみる.このサンプラーの特別な場合は既出の例で説明してある.ここでは,分布は2つのパラメータ m および α で定義される.

以下はギブスサンプラーを定義する.
In[114]:=
Click for copyable input

既出の特別なギブスサンプラーでは,m が16で α が2であった.

以下は のサンプラーからの5つのベクトルである.
In[115]:=
Click for copyable input
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): 721741.

[2] Casella, G. and E. I. George. "Explaining the Gibbs Sampler." The American Statistician 46, no. 3 (1992): 167174.

[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): 330.

[4] Nishimura, T. "Tables of 64-Bit Mersenne Twisters." ACM Transactions on Modeling and Computer Simulation 10, no. 4 (2000): 348357.

[5] Junod, P. "Cryptographic Secure Pseudo-Random Bits Generation: The BlumBlumShub 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): 341345.

[9] Cheng, R. C. H. and G. M. Feast. "Some Simple Gamma Variate Generators." Applied Statistics 28, no. 3 (1979): 290295.

[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): 515.

[12] Cheng, R. C. H. "Generating Beta Variables with Nonintegral Shape Parameters." Communications of the ACM 21, no. 4 (1978): 317322.

[13] Atkinson, A. C. "A Family of Switching Algorithms for the Computer Generation of Beta Random Variables." Biometrika 66, no. 1 (1979): 141145.

[14] Bailey, R. W. "Polar Generation of Random Variates with the t-Distribution." Mathematics of Computation 62, no. 206 (1994): 779781.

[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): 216223.

[17] Kachitvichyanukul, V. and B. W. Schmeiser. "Computer Generation of Hypergeometric Random Variates." Journal of Statistical Computation and Simulation 22, no. 2 (1985): 127145.

[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): 163179.

[19] Matsumoto, M. and T. Nishimura. "Dynamic Creation of Pseudorandom Number Generators", Monte Carlo and QuasiMonte Carlo Methods 1998, Springer, 2000, 5669.