Gamma分布からのサンプリング

自分用メモ。DPMにおいてハイパーパラメータのサンプリングをする必要がある場合、Gamma分布からサンプリングしてこなければならない。「Boostさんにお願いしてどうにかしてもらおう...」と思っていたところ、BoostさんはGamm分布の1変数のバージョンでしかAPIを提供していないということに気づいてしばらく止まる。

Gamma分布のpdf*1は、shape parameter k、scale parameter $\theta$ を用いて
$f(x; k, \theta) = x^{k-1} \frac{\exp{-\frac{x}{\theta}}}{\theta^k \Gamma(k)}$
と書ける。Boostさんはこのscale parameter $\theta$ が $\theta = 1$ の場合しか提供していないということである。しかし、shapeをあれこれ変えた場合のサンプリングも容易にできる。scaleのほうをfixすると
$f(x; k, \theta = 1) = x^{k-1} \frac{\exp{-x}}{\Gamma(k)}$
となるが、この場合で元の式で $t = \frac{x}{\theta}$ と変数変換すればこの式。ゆえにBoostさんで提供されているGamma分布でshapeを指定して乱数を生成し、その値をscaleで「かければ」両方指定したものと同じということである。まどろっこしい...。

プログラム

念のため実験して確かめる。Gamma分布の期待値は $\theta k$ であるが、それに近い感じになるかを確かめる。

#include <iostream>
#include <boost/random.hpp>

template<class D, class G = boost::mt19937>
class Rand {
  G gen_;
  D dst_;
  boost::variate_generator<G, D> rand_;
public:
  Rand() : gen_(static_cast<unsigned long>(time(0))), rand_(gen_, dst_) {}
  template<typename T1>
  Rand(T1 a1) : gen_(static_cast<unsigned long>(time(0))), dst_(a1), rand_(gen_, dst_) {}
  template<typename T1, typename T2>
  Rand(T1 a1, T2 a2) : gen_(static_cast<unsigned long>(time(0))), dst_(a1, a2), rand_(gen_, dst_) {}

  typename D::result_type operator()() { return rand_(); }
};

int main(int argc, char *argv[]) {
  double shape = 10.0;
  double scale = 2.0;

  Rand<boost::gamma_distribution<> > rgamma(shape);
  double result = 0.0;
  int N = 1000;
  for (int i = 0; i < N; i++) {
	result += rgamma() * scale;
  }
  std::cout << (double) result / N << std::endl;
}

それっぽい値になったので大丈夫そうである。

% g++ gamma.cpp; ./a.out
20.0174

注意

ひじょーにconfusingなんだけど、Gamma分布は
$f(x; \alpha, \beta) = x^{\alpha - 1} \frac{\beta^\alpha \exp^{- \beta x}}{\Gamma(\alpha)}$
とパラメタライズされることがある(PRMLなどはこっちの形式で書かれていて自分もこっちに慣れている)。この場合、上のコードで「かけているところ」は「割らないといけない」ことに注意しよう。というか、この違いではまって「変な乱数生成されてるけどなんで...」としばらくはまったのであった。。。