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

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

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

$$Y_1\sim Gamma({\alpha }_1,\theta ),...,Y_K\sim Gamma({\alpha }_K,\theta )$$

则有：

$$\begin{array}{rl}V& =\sum _{i=1}^K{Y}_i\sim Gamma({\alpha }_0,\theta ),\\ X& =(X_1,...,X_K)=(\frac{Y_1}{V},...,\frac{Y_K}{V})\sim Dir({\alpha }_1,...,{\alpha }_K)\end{array}$$

就是说，可以使用 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;
    }
}