質問:
離散分布からサンプリングする方法は?
Barry
2013-08-21 01:40:40 UTC
view on stackexchange narkive permalink

単一の確率変数Xからの可能な結果を​​支配する分布があると仮定します。これは、Xが1、2、3、4のいずれかの値である場合の[0.1、0.4、0.2、0.3]のようなものです。

この分布からサンプリングすることは可能ですか。つまり、その結果の確率が与えられた場合に、考えられる結果のそれぞれに疑似乱数を生成することは可能ですか。したがって、2を取得する確率を知りたい場合、サンプリング操作は0.34などを返す可能性があります。

私が尋ねる理由は、アクション選択ポリシーを実装しようとしているためです。研究論文に基づく強化学習法。私が論文から集めたものから、著者は「適応数値積分によって得られた累積確率密度関数を介して一様分布U [0,1]をマッピングする」ことによって分布をサンプリングすることができます。これから、彼は各試行の遷移確率をサンプリングします...

これに関する情報をいただければ幸いです...

よろしくお願いします

離散確率分布をサンプリングするには、さまざまな方法があります。このペーパーでは、累積分布関数を使用します($ u <0.1 $が「1」を出力する場合、$ <0.1 + 0.4 $が「2」を出力する場合など、(0,1)で$ U = u $を生成します)。速度が問題になる場合(たとえば、何十億回もサンプリングしたい場合)、はるかに効率的な方法があります。
@Glen_b離散rvをサンプリングするためのより効率的な方法を挙げてください。これはとても興味深いです。
@Rigaは以下の私の答えを参照してください
ここに「エイリアスメソッド」を説明する素晴らしい記事があります:http://www.keithschwarz.com/darts-dice-coins/
五 答え:
jtobin
2013-08-21 02:19:57 UTC
view on stackexchange narkive permalink

もちろんです。これは、その分布から n 回サンプリングし、置換するR関数です。

  sampleDist = function(n){sample(x = c(1,2、 3,4)、n、replace = T、prob = c(0.1、0.4、0.2、0.3))}#> sampleDist(10)#[1] 4 2 2 2 2 2 4 1 2 2  

もう少し低いレベルに移動したい場合は、Rソース(Cで記述)をチェックすることで実際に使用されているアルゴリズムを確認できます。

  / *不等確率サンプリング; with-replacement case * nはpとpermの長さです。 pには確率が含まれ、perm *には実際の結果が含まれ、ansにはサンプリングされた値の配列*が含まれます。 * / static void ProbSampleReplace(int n、double * p、int * perm、int nans、int * ans){double rU; int i、j; int nm1 = n-1; / *要素IDを記録します* / for(i = 0; i < n; i ++)perm [i] = i + 1; / *確率を降順で並べ替えます* / revsort(p、perm、n); / *累積確率を計算します* / for(i = 1; i < n; i ++)p [i] + = p [i-1]; / *サンプルを計算します* / for(i = 0; i < nans; i ++){rU = unif_rand(); for(j = 0; j < nm1; j ++){if(rU < = p [j])break; } ans [i] = perm [j]; }}  
さて、私は今何が起こっているのか理解しています。すべての答えに感謝します。これが他の誰かに役立つことを本当に望んでいます。私はすべての正しい答えを選択できればよかったのに。皆に感謝します
なぜ離散分布をソートする必要があるのですか?
@PavithranIyer私もそれを見ませんでした(または見ませんでした)。おそらくそれは最適化の試みです。最初に大きな確率をテストすることは、ループの早い段階でより頻繁に停止できることを意味します。しかし、頻繁にサンプリングしない限り(大きな「nans」)、最初のソートのコストに見合う価値があるとは思えません。しかし、繰り返しになりますが、数回サンプリングするだけでは、違いに気付くことはありません。
Glen_b
2013-08-22 08:17:49 UTC
view on stackexchange narkive permalink

コメントの質問に答えて、cdf法よりも離散分布を行う可能性のある*高速な方法の概要を以下に示します。

*離散ケースによってはうまくいくため、「潜在的に」と言います。実装された逆累積分布関数アプローチは非常に高速です。一般的なケースは、追加のトリックを導入せずに高速化するのは困難です。

質問の例のように4つの異なる結果の場合、逆累積分布関数アプローチ(または実質的に同等のアプローチ)のナイーブバージョンは次のようになります。罰金;ただし、カテゴリが数百(または数千、数百万)ある場合は、少し賢くならずに遅くなる可能性があります(最初のカテゴリが見つかるまで、累積分布関数を順番に検索したくないことは確かです。累積分布関数がランダムユニフォームを超える場合)。それよりも速いアプローチがいくつかあります。

[インデックスを使用して値を見つけるための順次よりも速いアプローチに関連している、以下で言及する最初のいくつかのことを見ることができます。単なる「cdfを使用したよりスマートなバージョン」です。もちろん、「ソートされたファイルの検索」などの関連する問題を解決するための「標準」アプローチを検討すると、シーケンシャルパフォーマンスよりもはるかに高速なメソッドになります。適切な関数を呼び出すことができれば、そのような標準的なアプローチで十分な場合があります。]

とにかく、離散分布から生成するためのいくつかの効率的なアプローチへ。


1 )「テーブルメソッド」。 $ n $ span>カテゴリの $ O(n)$ span>になる代わりに、一度設定すると、「シンプル"(a)のこれのバージョン(ディストリビューションが適切な場合)は $ O(1)$ span>です。


a)単純なアプローチ-有理確率を仮定する(上記のデータ例で実行):
-「1」、4つの「2」、2つの「3」、および3つの「4」を含む10個のセルで配列を設定します。離散一様分布(連続一様分布から簡単に実行可能)を使用してサンプルを作成すると、単純で高速なコードが得られます。

b)より複雑なケース-「良い」確率は必要ありません。 $ 2 ^ k $ span>セルを使用します。そうしないと、使用するセルがそれよりも少なくなります。したがって、たとえば、次のことを考慮してください。

  x 0 1 2 3 4 5 6P(X = x)0.4581 0.0032 0.1985 0.3298 0.0022 0.0080 0.0002  

(もちろん、10000個のセルを持ち、以前の正確なアプローチを使用することもできますが、これらの確率が不合理である場合はどうでしょうか?)

$ k = 8 $を使用しましょう span>。確率に $ 2 ^ k $ span>を掛け、切り捨てて、必要な各セルタイプの数を調べます。

  x 0 1 2 3 4 5 6 TotP(X = x)0.4581 0.0032 0.1985 0.3298 0.0022 0.0080 0.0002 1.0000 [256p(x)] 117 0 50 84 0 2 0 253  

次に、最後の3つのセルは基本的に「代わりに、この他の分布から生成する」(つまり、p(x)-\ frac {\ lfloor 256 p(x)\ rfloor} {256}をpmfに正規化):

  x * 0 1 2 3 4 5 6 TotP(X = x *)0.091200 0.273067 0.272000 0.142933 0.187733 0.016000 0.017067 1.000000  

「波及効果」テーブルは、任意の合理的な方法で実行できます(ここに到達する時間は、約1%であり、それほど高速である必要はありません)。したがって、 $ \ frac {253} {256} $ span>の時間のランダムユニフォームを生成し、最初の8ビットを使用してランダムセルを選択し、次の値を出力します。セル;初期設定後、これらすべてを非常に高速に行うことができます。 「2番目のテーブルから生成」というセルにヒットしたときのもう1つの $ \ frac {3} {256} $ span>。ほとんどの場合、 $(0,1)$ span>で単一のユニフォームを生成し、乗算、切り捨て、および配列要素へのアクセスのコストから離散乱数を取得します。

2)「ヒストグラムの2乗」法。これは(1)に関連していますが、各セルは、(連続)一様に応じて、実際には 2つの値のいずれかを生成できます。したがって、1からnまでの離散値を生成し、それぞれの中で、そのメイン値と2番目の値のどちらを生成するかを確認します。有界確率変数で動作します。波及効果テーブルはなく、通常、方法(1)よりもはるかに小さいテーブルを使用します。通常、1:nを選択すると、一様乱数の最初の数ビットが使用されるように設定され、残りの部分は、そのビンの2つの値のどちらを出力するかを示します。

メソッドの概要を説明する最も簡単な方法は、上記の例でそれを行うことです。

分布を4つのビンを持つヒストグラムと考えてください:

original 'histogram'

最も高いバーの上部を切り取り、短いバーに配置して、「二乗」します。バーの平均「高さ」は0.25になります。したがって、2番目のバーを0.15カットして最初のバーに配置し、0.05を4番目のバーからカットして3番目のバーに配置します。

'squaring off' the histogram

これを1つの色が複数のビンになる可能性がありますが、ビンが2色を超えることはありません。

ここで、4つのビンの1つをランダムに選択します(ユニフォームの上から2つのランダムビットが必要です)。次に、残りのビットを使用して、均一に分散された垂直位置を指定し、色間の区切りと比較して、2つの値のどちらを出力するかを決定します。非常に高速ですが、通常は「テーブル」メソッドほど高速ではありません。

-

これらのメソッドは、無制限の変数を処理するように適合させることができますが、ここでも、ほとんどの場合高速です。 '。

参照: http://www.jstatsoft.org/v11/i03/paper

これらの比較的遅い部分は、値の表;何を生成するかがわかっている場合(「将来、この分布から値を何度もサンプリングする必要がある」)、作成しようとするよりも適しています。 「このASAPから100万の値をサンプリングする必要がありますが、再度実行する必要はありません」と、さまざまな優先順位が作成されます。多くの場合、ソートされた値を検索するための「標準的なコンピューティングアプローチ」(つまり、cdfメソッドをより迅速に実行するため)のいくつかは、実際には最善の選択である可能性があります。


他にも高速なアプローチがあります離散分布からの生成に。注意深くコーディングすれば、非常に高速な生成を行うことができます。例:

3)棄却法(「棄却-棄却」)は、離散分布を使用して実行できます。すでに高速な方法で生成できるスケールアップされた離散pmfである離散メジャー化関数(「エンベロープ」)がある場合、それは直接適応し、場合によっては非常に高速になる可能性があります。より一般的には、連続分布から生成できることを利用できます(たとえば、結果を離散化して離散エンベロープに戻すことにより)。

ここで、便利な累積分布関数(または逆累積分布関数)がない離散確率関数 $ f $ span>があると想像してください- -実際、この図では正規化定数すら持っていなかったため、プロットは正規化されていません:

plot of an unnormalized unimodal discrete probability mass function on the natural numbers; it has its mode at 2 and eventually tails off approximately geometrically

次に、離散確率関数 $ g $ span>から生成するのに便利な関数を見つける必要があります。これには定数 $ c $ span>であり、少なくとも $ f $ span>と同じ大きさである必要があります(これがすべての $ x $ span>値)。つまり、すべての可能な $ x $ span>に対して $ cg(x)\ geq f(x)$ span>です。値。

適切な $ g $ span>を簡単に特定できる場合もありますが、1つの便利なオプションは、左側の部分に離散一様分布を混合することです。右側の $ f $ span>と少なくとも同じくらい重いテールを持つ分布。そのための2つの合理的に便利な選択肢は、幾何分布(テールが指数関数よりもゆっくりと減少していない場合)と、 $ \を取ることによって得られる離散化パレートまたは離散化ハーフコーシー分布のようなものです。一部のパレートまたは半コーシーランダム変量 $ X $ span>のlfloorX \ rfloor $ span>(いずれの場合も、pmfが指数関数的よりもゆっくりと減少している場合)。

(さらに言えば、幾何学自体は指数を離散化することで生成できます。)

この場合、左側の離散一様分布と右側の幾何学が非常にうまく機能します。 :

The previous discrete pmf with the aforementioned uniform&geometric envelope (majorizing function)

(注意:ここにプロットされているのは正規化されていないpmfであるため、y軸は確率を表すのではなく、比例するものを表します。確率に)

次に、手順は、 $ g $ span>から提案された値 $ x $ span>をシミュレートすることです。 $(0、cg(x))$ span>で $ U $ span>をシミュレートし、 $ U<f $ span>、提案された $ x $ span>を受け入れます(それ以外の場合は拒否し、新しく提案された $ x $ span>)。

ありがとう、グレン! $ 2 ^ k $アプローチは有望です。例を挙げて、「ヒストグラムの二乗」方法を明確にしていただけませんか。
確かに、私は行くことができますが、あなたの元のコメントが私に*名前*アプローチを求めたことを思い出してください...そしてそれがその名前です。それまでの間、基本的なプロセスについての適切な説明があります[ここ](http://www.robertowor.com/csci4151/lecture3.htm)。 [ロビンフッド法](http://www.jstatsoft.org/v11/i03/paper)とも呼ばれます。
@Riga質問の例の2番目のケースのアイデアの簡単な概要で説明を更新しました。
グレン、ありがとうございました!あなたの参照と説明は貴重な知識です。
こんにちは@Glen_b,は、2番目/波及効果のテーブル確率をどのように思いついたのか詳しく教えてください。
@greendiod $ p(x)-\ frac {\ lfloor 256 p(x)\ rfloor} {256} $は合計が1にスケールアップされるため、この例では、256/3が乗算されます。
@Glen_bわかりました、なるほど、あなたは不足している確率を埋め、[0,1(。ランダムユニフォームの8ビット?
確率変数の迅速な生成を実装するのに適したほとんどの言語は、ビット演算を提供します。したがって、Cで作業している場合(または、Cコンパイラやその他のコンパイラが不十分なときに何度も行った、チップセットを低レベルにするためのアセンブリ言語での記述)は、次のようになります。標準。何らかの理由で高速ビット操作が不足している非常に高水準の言語で書き込もうとしている場合は、おそらくその程度までパフォーマンスを最適化することにあまり関心がないでしょう。別のユニフォームを使用するだけです。... ctd
ctd ...一方、「*この*言語でビット操作を行うにはどうすればよいですか」のような質問をしている場合、その質問には間違ったフォーラムです。
@Glen_b実際、私の質問は概念レベルにあります(いつものように、詳細には悪魔がいます)。ビットレベルの演算子がなかった場合は、$ U([0、...、255])$に続く古典的な$ floor(U * 256)$によってランダムなセルを選択していました(あなたのように)遅い場合でも、「波及効果」テーブルから生成することを提案しました)。さて、$ U $は、Cで言うと、IEEE754ダブル(実装の詳細は大丈夫です)で、$ U $の最初の8ビットで十分ですか?(なぜ?リンクは?)
@green均一なRNGは通常、実際には浮動小数点数ではなく$ 0,1、...、m-1 $(一部の$ m $の場合)の整数を生成し、それらを$の浮動小数点形式(doubleなど)にスケーリングします。[0,1)$($ m $で割ることによる)。doubleの小数部分の上位ビット(または他の浮動小数点形式が使用されているもの)ではなく、整数の上位ビットを使用するのが最も簡単です。一部のRNGでは十分にランダムではないため、ローエンドビットは避けます。
@Glen_b波及効果テーブルから提案された世代に惑わされました。実際、Rnd(0、m -1)からU、Rnd(0、n -1)に移動する必要がありました。すべての入力に感謝します
通常、$ m $は$ n $と比較して巨大になります
$ U
@Fooはい。あいまいまたは明確さが不十分だと思われる場合は、編集できます。
@Glen_bプロットでは、縦軸のラベルが誤っている必要があります。離散確率変数は、単一の観測で1より大きいメジャーを持つことはできません。すべての値の合計は1になる必要があり、それらはすべて負ではありません。
@Lucas誤ったラベルは付けられていません-私の答えは、プロットされるものが*非正規化* $ f $であると明示的に述べています。ここで、「*この図では、正規化定数すら持っていなかったので、プロットは非正規化されています*」(つまり、変数pxは実際の確率に比例するだけです)。変数 `px`を確率を表すと解釈すると、私の答えを誤解することになります。変数pxがpmf($ f $)の正規化されていないバージョンを表すという説明を、私が推測する答えのどこかでもう一度繰り返すことができます。
@Glen_bあなたは正しいです、答えは密度がスケールアップされていると述べています-私は誤解しています-答えが長いとすれば、答えのどこか下に言い換えると理解しやすくなるかもしれません。誤読をお詫びします。
ありがとう;私は今そうしました-最後のプロット(メジャー化関数が含まれているもの)の下に文を挿入しました。( `qx`のようなものではなく` px`を使用したので)変数名の選択は、これが実際の確率であるという印象に貢献した可能性があります。
user1893354
2013-08-21 01:51:40 UTC
view on stackexchange narkive permalink

Pythonでは、scipy.statsから

 のようなことを行うことができますimportrv_discrete x = [1,2,3,4] px = [0.1,0.4,0.2,0.3] sample = rv_discrete(values =(x、px))。rvs(size = 10) 

これにより、分布から10個のサンプルが得られます。これを繰り返して、2であるサンプルの比率を見つけることができます。

Greg Snow
2013-08-21 02:22:14 UTC
view on stackexchange narkive permalink

はい、それは可能でかなり簡単です。正確には、使用しているツールによって異なります。

Rでは、 sample(1:4、n、prob = c)になります。 (0.1,0.4,0.2,0.3)、replace = TRUE)ここで、 n はサンプリングする値の数です。

同等の機能のないツールの場合それでも均一な値を生成でき、RVは0.1未満の場合は1、0.1から0.5の場合は2、0.5から0.7の場合は3、0.7より大きい場合は4になります(これは、にマッピングするという考え方です。累積)。

この例では、1、4、2、2、3、3の4のセットから均一にサンプリングして、同じ確率を取得することもできます。

+1。 Excelの同等の関数の実用的な例として(これがいかに簡単かを示すために)、確率の累積合計の配列を設定し(たとえば、(0,0.1,0.5,0.7,1.0))、名前を付けます(例: `cumsum`)、` = MATCH(RAND()、cumsum、1) `を使用して、確率0.1、0.4、0.2、0.3の値1、2、3、4を生成します。これは、重み付けされたサンプリングが配列検索にどのように関連しているかを明確に示しています。
この答えを本当にありがとう、それは本当に私を大いに助けました、私はこれが他の誰かを助けることを本当に望んでいます。私はすべての正しい答えを選択できればよかったのに。皆に感謝します
Nick Cox
2013-08-21 02:40:41 UTC
view on stackexchange narkive permalink

Stataの場合:

Mataの場合、 http://www.stata.com/help.cgi?mf_runiform rdiscrete()を使用します。 >

Stata自体には、さまざまな方法があります。これが1つです:

 。 gen rnd = runiform()。 gen y = cond(rnd < = 0.1、1、cond(rnd < = .5、2、cond(rnd < = .7、3、4))) 


このQ&Aは英語から自動的に翻訳されました。オリジナルのコンテンツはstackexchangeで入手できます。これは、配布されているcc by-sa 3.0ライセンスに感謝します。
Loading...