確率統計
超幾何分布サンプリング; Hyper Geometric Distribution Sampling

概要

---

超幾何分布サンプリング
超幾何分布とは以下の式で表現される確率分布。ある属性をもつ\(K\)個の要素を含む\(N\)の要素からなる母集団から\(M\)個の要素を抽出した時、 その属性をもつ要素の個数が\(k\)である確率を示す。なお\(K\)と\(M\)を入れ替えても同じ確率分布になる。

$$ P[X=k] = \frac{ {}_K \mathrm{C}_k \dot {}_{N-K} \mathrm{C}_{M-k} }{ {}_N \mathrm{C}_M }, \quad Max(0, M+K-N) \leq k \leq Min(M, K) $$


超幾何分布に従う確率変数の生成にはルーレット選択アルゴリズムを用いる。そのために確率分布を求めるが、二項係数から求めようとするとオーバーフローなどの問題が生じる。そこで以下の漸化式を用いる。
$$ \begin{eqnarray} p(k; N, K, M) &=& \frac{ {}_K \mathrm{C}_k \dot {}_{N-K} \mathrm{C}_{M-k} }{ {}_N \mathrm{C}_M } \\ p(M; N, N, M) &=& 1 \\ p(k; N, K-1, M) &=& \frac{K-k}{K} p(k; N, K, M)+\frac{k+1}{K}p(k+1; N, K, M) \end{eqnarray} $$
ソースコード

---

namespace ExRandom.Discrete {
    public class HyperGeometricRandom : Random {
        readonly MT19937 mt;
        readonly int[] indexs;
        readonly double[] probs;

        public HyperGeometricRandom(MT19937 mt, int n = 30, int k = 20, int m = 10) { 
            if(mt == null) {
                throw new ArgumentNullException();
            }

            if(n <= 0 || k <= 0 || m <= 0 || n < k || n < m) {
                throw new ArgumentException();
            }

            if(k < m) {
                int swap;
                swap = k;  k = m;  m = swap;
            }

            int max = m, min = Math.Max(0, k + m - n);
            double p;
            double[] next_probs;

            this.mt = mt;

            this.indexs = new int[max - min + 1];
            this.probs = new double[max - min + 1];

            for(int i = 0; i < indexs.Length; i++) {
                indexs[i] = i + min;
            }

            probs[max - min] = 1;

            for(int i = n, j; i > k; i--) {
                next_probs = new double[max - min + 1];

                for(j = Math.Max(max - n + i, min); j <= max; j++) {
                    if(j > 0) {
                        p = (double)(i - j) / (double)(i);

                        next_probs[j - min] = probs[j - min] * p;
                        next_probs[j - min - 1] += probs[j - min] * (1 - p);
                    }
                    else { 
                        next_probs[j - min] = probs[j - min];
                    }
                }

                probs = next_probs;
            }

            Array.Sort(this.probs, this.indexs);
        }

        public override int Next() {
            double r = mt.NextDouble_OpenInterval1();

            for(int i = probs.Length - 1; i >= 0; i--) {
                r -= probs[i];
                if(r <= 0) {
                    return indexs[i];
                }
            }

            return indexs[probs.Length - 1];
        }
    }
}

関連項目

---

メルセンヌ・ツイスタ
各種確率分布サンプリング基本クラス