C++でディリクレ分布からサンプルする

はじめに

最近「続・わかりやすいパターン認識」という本を読んでいます。いわゆるノンパラベイズ本でして、パラメータの事前分布としてディリクレ分布が多用されています。気になったのでC++でサンプリングする方法を調べてみました。

ディリクレ分布のサンプリング

Googleに"dirichlet sampling c++"などと入れてライブラリを探してみると、GNU Scientific Library (GSL)がおすすめと書いてありました。GPLです。
random - Looking for C/C++ implementations of sampling from multinomial and Dirichlet distributions - Computational Science Stack Exchange

他にライブラリを探しても良いのが見つからなかったので困っていましたが、ガンマ分布のサンプルを合計が１になるように正則化すると、ディリクレ分布からのサンプリングと等価となるようです。なぜかは書いてないので分かりませんが、簡単そうなので実装してみます。

sampling - Drawing from Dirichlet distribution - Cross Validated

ガンマ分布

ガンマ分布の実装はC++11標準で対応しています。最初の引数でalphaを指定できます。

std::default_random_engine engine;
std::gamma_distribution<double> gamma(7.5, 1.0);

std::ofstream ofs("result.csv");
for (int i = 0; i < 1000; i++)
{
  ofs << gamma(engine) << std::endl;
}
ofs.close();

ディリクレ分布

ガンマ分布の値を合計して割ることで生成できます。
三次元で表示するのがむずかしかったので、サンプルした最初の二要素をプロットしています。C++に使いやすい数値計算ライブラリがないのは何でなんでしょうね。SciPyとかR使えば良いんですけど。

std::default_random_engine engine; // 同じseed
std::gamma_distribution<double> gamma;
VectorXd DirichletSampling(const VectorXd& alpha)
{
  VectorXd x(alpha.size());
  for (int i = 0; i < alpha.size(); i++)
  {
    gamma.param(alpha(i));
    x(i) = gamma(engine);
  }
  x /= x.sum();

  return x;
}

int main()
{
  VectorXd alpha(3);
  alpha << 2.0, 4.0, 6.0;

  std::ofstream ofs("result.csv");
  for (int i = 0; i < 10000; i++)
  {
    VectorXd x = DirichletSampling(alpha);
    ofs << x(0) << "," << x(1) << "," << x(2) << std::endl;
  }
  ofs.close();
}