ディリクレ分布サンプリング - わかばテクノロジ

確率統計
ディリクレ分布サンプリング; Dirichlet Distribution Sampling

概要

$K-1$次ディリクレ分布は以下の式で表現される確率分布。ここで$\boldsymbol{\alpha}$は形状パラメータである。ベータ分布を多変量に拡張した分布である。

$$f(\boldsymbol{x})=\frac{\Gamma \left(\sum_{i=1}^{K} \alpha_i \right)}{\prod_{i=1}^{K} \Gamma(\alpha_i)} \displaystyle \prod_{i=1}^{K} x_i^{\alpha_i - 1} $$

ディリクレ分布に従う確率変数の生成はベータ分布と同様に、形状パラメータ$\kappa=\alpha_i $、尺度パラメータ$\theta=1$の$K$個のガンマ分布に従う確率変数を$R_i$としたとき、$\boldsymbol{x} = (\frac{R_1}{\sum R_i}, \ldots, \frac{R_{K-1}}{\sum R_i})^T$がディリクレ分布に従うことを利用する。

ソースコード


namespace ExRandom.MultiVariate {
    public class DirichletRandom : Random<double> {
        readonly int dim;
        readonly Continuous.GammaRandom[] gr_list;

        public DirichletRandom(MT19937 mt, params double[] alpha_list) {
            if(mt == null) {
                throw new ArgumentNullException();
            }

            if(alpha_list == null || alpha_list.Length <= 1) {
                throw new ArgumentException();
            }

            this.dim = alpha_list.Length;

            this.gr_list = new Continuous.GammaRandom[dim];

            for(int i = 0; i < dim; i++) {
                this.gr_list[i] = new Continuous.GammaRandom(mt, kappa:alpha_list[i], theta:1);
            }
        }

        public override Vector<double> Next() {
            double r_sum = 0;
            double[] r_list = new double[dim];
            Vector<double> v = new Vector<double>(new double[dim - 1]);

            for(int i = 0; i < dim; i++) {
                r_list[i] = gr_list[i].Next();
                r_sum += r_list[i];
            }

            r_sum = Math.Max(r_sum, double.Epsilon);

            for(int i = 0; i < dim - 1; i++) {
                v[i] = r_list[i] / r_sum;
            }

            return v;
        }
    }
}