質問

シミュレーションシステムに取り組んでいます。すぐに、いくつかのシミュレーション入力の値の実際の分布の実験データ(ヒストグラム)を取得します。

シミュレーションの実行時に、測定された分布に一致するランダムな値を生成できるようにしたいと思います。元のヒストグラムを保存せずにこれを実行したいです。いくつかの良い方法は何ですか

  1. ヒストグラムを分布を表す一連のパラメーターにマッピングしますか?
  2. 実行時にこれらのパラメーターに基づいて値を生成しますか?

編集:入力データは、いくつかの異なるタイプのイベントの期間です。異なるタイプには異なる分布関数があると予想されます。

役に立ちましたか?

解決

少なくとも2つのオプション:

  1. ヒストグラムを統合し、数値的に反転します。
  2. 拒否

数値統合

ウィリアム・R・ギブスによる現代物理学の計算より:

  

常に[関数]を数値的に積分し、[ cdf ]を反転させることができます   しかし、特に pdf が変化している場合、これはあまり満足のいくものではありません。   急速に。

文字通り、範囲 [0-1)をターゲット分布の適切な範囲に変換するテーブルを構築します。次に、通常の(高品質の)PRNGをスローして、テーブルで翻訳します。面倒ですが、明確で実行可能で、完全に一般的です。

拒否:

ターゲットヒストグラムを正規化してから、

  1. サイコロを投げて、範囲に沿ってランダムに位置( x )を選択します。
  2. もう一度投げて、新しい乱数がこのビンの正規化されたヒストグラムよりも小さい場合は、このポイントを選択します。それ以外の場合は(1)に進みます。

繰り返しますが、心はシンプルですが、明確で機能しています。非常に低い確率(長い尾を持つピーク)での配布には時間がかかる場合があります。


これらの方法の両方で、ステップ関数ヒストグラムが望ましくない場合、区分多項式近似またはスプラインでデータを近似して滑らかな曲線を生成することができます---後で残すそれは時期尚早な最適化かもしれません。


特殊なケースには、より良いメソッドが存在する場合があります。

これらはすべて非常に標準的なものであり、詳細が必要な場合は数値分析の教科書に表示されるはずです。

他のヒント

問題に関する詳細情報が役立ちます。たとえば、ヒストグラムはどのような値になりますか?カテゴリ(色、文字など)か、連続(高さ、時間など)ですか?

ヒストグラムがカテゴリデータを超えている場合、カテゴリ間に多くの相関関係がない限り、分布をパラメータ化することは難しいと思います。

ヒストグラムが連続データを超える場合、ガウス分布の混合を使用して分布を近似しようとする場合があります。つまり、$ \ sum_ {i = 1} ^ n w_i N(m_i、v_i)$を使用してヒストグラムを近似しようとします。ここで、m_iおよびv_iは平均と分散です。次に、データを生成する場合は、まず重みw_iに比例する確率で1..nからiをサンプリングし、次に任意のガウスから行うようにx〜n(m_i、v_i)をサンプリングします。

どちらの方法でも、混合モデルの詳細を読むことができます。

そのため、特定の確率分布を生成するために必要なのは、分位数関数、これはの逆です @dmckeeが言うように、累積分布関数

質問は次のようになります:与えられた連続ヒストグラムを記述する変位値関数を生成して保存する最良の方法は何ですか?答えは入力の形状に大きく依存すると感じています-それが何らかのパターンに従う場合、最も一般的なケースを単純化する必要があります。ここに行くときに更新します。


編集:

今週、この問題を思い出させる会話がありました。ヒストグラムを方程式として記述するのを忘れて、テーブルを保存するだけであれば、O(1)時間で選択を行うことができますか?精度を損なうことなく、O(N lgN)の構築時間を犠牲にしてできることがわかりました。

N個のアイテムの配列を作成します。配列への一様なランダム選択は、確率が1 / Nのアイテムを見つけます。各項目について、この項目が実際に選択されるヒットの割合と、この項目が選択されない場合に選択される別の項目のインデックスを格納します。

重み付きランダムサンプリング、C実装:

//data structure
typedef struct wrs_data {
  double share; 
  int pair;
  int idx;
} wrs_t;


//sort helper
int wrs_sharecmp(const void* a, const void* b) {
  double delta = ((wrs_t*)a)->share - ((wrs_t*)b)->share;
  return (delta<0) ? -1 : (delta>0);
}


//Initialize the data structure
wrs_t* wrs_create(int* weights, size_t N) {
  wrs_t* data = malloc(sizeof(wrs_t));
  double sum = 0;
  int i;
  for (i=0;i<N;i++) { sum+=weights[i]; }
  for (i=0;i<N;i++) {
    //what percent of the ideal distribution is in this bucket?
    data[i].share = weights[i]/(sum/N); 
    data[i].pair = N;
    data[i].idx = i;
  }
  //sort ascending by size
  qsort(data,N, sizeof(wrs_t),wrs_sharecmp);

  int j=N-1; //the biggest bucket
  for (i=0;i<j;i++) {
    int check = i;
    double excess = 1.0 - data[check].share;
    while (excess>0 && i<j) {
      //If this bucket has less samples than a flat distribution,
      //it will be hit more frequently than it should be.  
      //So send excess hits to a bucket which has too many samples.
      data[check].pair=j; 
      // Account for the fact that the paired bucket will be hit more often,
      data[j].share -= excess;  
      excess = 1.0 - data[j].share;
      // If paired bucket now has excess hits, send to new largest bucket at j-1
      if (excess >= 0) { check=j--;} 
    }
  }
  return data;
}


int wrs_pick(wrs_t* collection, size_t N)
//O(1) weighted random sampling (after preparing the collection).
//Randomly select a bucket, and a percentage.
//If the percentage is greater than that bucket's share of hits, 
// use it's paired bucket.
{
  int idx = rand_in_range(0,N);
  double pct = rand_percent();
  if (pct > collection[idx].share) { idx = collection[idx].pair; }
  return collection[idx].idx;
} 

編集2:   少し調査した結果、O(N)時間で構築を行うことさえ可能であることがわかりました。注意深い追跡により、大小のビンを見つけるために配列をソートする必要はありません。 更新された実装はこちら

離散ポイントの重み付き分布で多数のサンプルを取得する必要がある場合は、同様の質問への回答

ただし、ヒストグラムを使用して連続したランダム関数を近似する必要がある場合、おそらくdmckeeの数値積分の答えが最善の策です。または、エイリアスを使用して、ポイントを左に保存し、2つのポイントの間で統一された数を選択することもできます。

ヒストグラムから選択するには(オリジナルまたは縮小)、 Walkerのエイリアスメソッド 高速かつシンプルです。

正規分布の場合、以下が役立つ場合があります。

http://en.wikipedia.org/wiki/Normal_distribution#Generating_values_for_normal_random_variables

ライセンス: CC-BY-SA帰属
所属していません StackOverflow
scroll top