C++ 迪利克雷（Dirichlet）分布

遇到一个要使用 dirichlet 分布的情形，发现 C++ 标准库中没有现成的。查阅维基百科发现，虽然它挺复杂，但是它跟 Gamma 分布有如下关系：

设有 K 个相互独立且分别满足 Gamma 分布的分布：

\[Y_1 \sim Gamma(\alpha_1, \theta), ..., Y_K \sim Gamma(\alpha_K, \theta) \]

则有：

\[\begin {aligned} V &= \sum_{i=1}^{K}Y_i \sim Gamma(\alpha_0, \theta),\\ X &= (X_1, ..., X_K) = \left(\frac {Y_1}{V}, ..., \frac {Y_K}{V}\right) \sim Dir(\alpha_1, ..., \alpha_K) \end {aligned} \]

就是说，可以使用 gamma 分布来生成 dirichlet 分布。以下是用 C++ 实现的版本：

#include <random>

void dirichlet(double* out, int k, double a) {
    std::gamma_distribution<double> gamma(a);
    std::random_device rd;
    double *y = alloca(sizeof(double)*k); // stack allocation
    double sum=0;
    for (int i=0; i<k; ++i) {
        y[i] = gamma(rd);
        sum += y[i];
    }
    for (int i=; i<k; ++i) {
        out[i] = y[i]/sum;
    }
}

void dirichlet(double* out, int k, double const* a) {
    using Gamma = std::gamma_distribution<double>;
    Gamma gamma;
    std::random_device rd;
    double *y = alloca(sizeof(double)*k); // stack allocation
    double sum=0;
    for (int i=0; i<k; ++i) {
        y[i] = gamma(rd, Gamma::param_type(a[i], 1));
        sum += y[i];
    }
    for (int i=0; i<k; ++i) {
        out[i] = y[i]/sum;
    }
}