確率統計
二項分布サンプリング; Binomial Distribution Sampling

概要

---

二項分布は以下の式で表現される。成功確率\(p\)である独立な試行を\(n\)回行ったとき\(k\)回成功する確率を示したものである。

$$ P[X=k] = \binom{n}{k} p^k (1-p)^{n-k} , \quad 0 \leq k \leq n $$


二項分布に従う確率変数の生成には、ルーレット選択アルゴリズムを用いる。これにより計算コストが比較的大きい乱数生成を1回にすることができる。
なお、ソースコード上では\(n\)が1024より大きいとき、計算誤差を避けるためルーレット選択を分割して行っている。

ソースコード

---

namespace ExRandom.Discrete {
    public class BinomialRandom : Random {
        const int skips = 1024;

        readonly RouletteRandom rd, rd_skip;
        readonly int n;

        public BinomialRandom(MT19937 mt, double prob = 0.5, int n = 10) {
            if(mt == null) {
                throw new ArgumentNullException();
            }

            if(n <= 0) {
                throw new ArgumentException();
            }
            
            this.n = n;

            if(n <= skips) {
                double[] p = Binomial.Coef(n, prob);
                this.rd = new RouletteRandom(mt, p);
            }
            else { 
                double[] p = Binomial.Coef(n % skips, prob);
                this.rd = new RouletteRandom(mt, p);

                double[] p_skip = Binomial.Coef(skips, prob);
                this.rd_skip = new RouletteRandom(mt, p_skip);
            }
        }

        public BinomialRandom(MT19937 mt, decimal prob, int n = 10) : this(mt, (double)prob, n){
        }

        public override int Next() {
            int i = n, cnt = 0;

            while(i > skips) {
                cnt += rd_skip.Next();
                i -= skips;
            }

            cnt += rd.Next();

            return cnt;
        }
    }
}

namespace ExRandom {
    internal static class Binomial {
        public static double[] Coef(int n, double p) {
            double[] b0 = {1}, b1;

            for(int i = 1, j; i <= n; i++) {
                b1 = new double[i + 1];

                for(j = 0; j < i; j++) {
                    b1[j] += (1 - p) * b0[j];
                    b1[j + 1] += p * b0[j]; 
                }

                b0 = b1;
            }

            return b0;
        }
    }
}

関連項目

---

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