サンプリングアルゴリズム

statistics

伝承サンプリング

β分布

ベータ分布 $\mathcal{\beta}(\alpha,\beta)$ からのサンプリング

$\alpha$ と $\beta$ が共に自然数の場合

std::mt19937 engine(std::random_device("")());
std::uniform_real_distribution<T> runif;

std::vector<T> u(alpha + beta - 1);
for (std::size_t i = 0, size = u.size(); i < size; ++i) {
  u[i] = runif(engine);
}
std::sort(u.begin(), u.end());

return u[alpha];

それ以外の場合

std::mt19937 engine(std::random_device("")());
std::gamma_distribution<T> rgamma1(alpha, T(1));
std::gamma_distribution<T> rgamma2(beta , T(1));

theta1 = rgamma1(engine);
theta2 = rgamma2(engine);

return theta1 / (theta1 + theta2);

データ拡大サンプリング

データ拡大サンプリング (data augmentation sampling)

あらまし

ある分布 $f(\vec{x})$ からサンプルを生成したいが、直接は難しい。 そこで、ダミー変数 $\vec{y}$ を用意して

f(\vec{x}) \propto \int_{\Omega_{\vec{y}}} f(\vec{x}, \vec{y}) d\mu(\vec{y})

となる $f(\vec{x}, \vec{y})$ をつくる。ここで、$f(\vec{x}, \vec{y})$ は $f(\vec{x} \mid \vec{y})$ および $f(\vec{y}, \vec{x})$ からのサンプリングが容易なように設計する。また、$\Omega_{\vec{y}}$ は $\vec{y}$ の定義域である。

アルゴリズム

$f(\vec{x}, \vec{y})$ からギブスサンプリングすれば所望のサンプルが得られる。すなわち、

  1. 初期サンプル $(\vec{x}, \vec{y})$ を設定
  2. $f(\vec{x} \mid \vec{y})$ から $\vec{x}$ をサンプリング
  3. $f(\vec{y} \mid \vec{x})$ から $\vec{y}$ をサンプリング

とする。系列

[\vec{x}_{i}, \vec{y}_{i}]_{i=1}^{N}

[\vec{x}_{i}]_{i=1}^{N}

部分だけに着目すると、$f(\vec{x})$ に従うサンプルになっている。

スライスサンプリング

あらまし

スライスサンプリング (slice sampling) は

f(x, y) = 1_{f(x) > y}

としたデータ拡大サンプリングである。ここで $1_{z}$ は $z$ が真のとき $1$ になり、偽のとき $0$ となる関数である。

\begin{aligned}
\int_{\Omega_{y}} f(x,y) d\mu(y)
  &= \int_{0}^{f(x)} dy \\
  &= f(x)
\end{aligned}

となるので条件を満たしている。

$f(x,y)$ の形を考えると理解しやすい。x-y 平面上でグラフに囲まれている図形を一様サンプリングするイメージ。

アルゴリズム

  1. 初期値を設定
  2. $y$ を $[0, f(x)]$ から一様サンプリング
  3. $f(x) > y$ となる区間から $x$ を一様サンプリング

連続の場合、基本的に単峰性の関数にしか使えない。

[0, 1]

$f(x)$ をサンプリングする。

double t = 0.5;
double x = t;
double y = f(x);
for (int i = 0; i < N; ++i) {
  // sample a slice
  std::uniform_real_distribution<> rslice(0, y);
  double const slice = rslice(engine);

  // chose new sample
  double left = 0.0;
  double right = 1.0;
  for (;;) {
    std::uniform_real_distribution<> runif(left, right);
    double const new_t = runif(engine);
    double const new_x = new_t;
    double const new_y = f(new_x);

    if ((new_y > slice) || (t == new_t)) {
      t = new_t;
      x = new_x;
      y = new_y;
      break;
    }

    if (t < new_t) {
      right = new_t;
    } else {
      left = new_t;
    }
  }

  // record x
  std::cout << x << '\n';
}

[0, ∞]

$f(x)$ をサンプリングする。$x = \frac{t}{1-t}$ という変数変換を使う。ヤコビアンは $dx = \frac{1}{(1-t)^{2}} dt$ となる。

double t = 0.5;
double x = t / (1 - t);
double y = f(x) / ((1 - t) * (1 - t));
for (int i = 0; i < N; ++i) {
  // sample a slice
  std::uniform_real_distribution<> rslice(0, y);
  double const slice = rslice(engine);

  // chose new sample
  double left = 0.0;
  double right = 1.0;
  for (;;) {
    std::uniform_real_distribution<> runif(left, right);
    double const new_t = runif(engine);
    double const new_x = new_t / (1 - new_t);
    double const new_y = f(new_x) / ((1 - new_t) * (1 - new_t));

    if ((new_y > slice) || (t == new_t)) {
      t = new_t;
      x = new_x;
      y = new_y;
      break;
    }

    if (t < new_t) {
      right = new_t;
    } else {
      left = new_t;
    }
  }

  // record x
  std::cout << x << '\n';
}

[-∞, ∞]

$f(x)$ をサンプリングする。$x = -\log\left(\frac{1}{t} - 1\right)$ という変数変換を使う。ヤコビアンは $dx = \frac{1}{t (1-t)} dt$ となる。

double t = 0.5;
double x = std::log(1.0 / t - 1.0);
double y = f(x) / (t * (1 - t));
for (int i = 0; i < N; ++i) {
  // sample a slice
  std::uniform_real_distribution<> rslice(0, y);
  double const slice = rslice(engine);

  // chose new sample
  double left = 0.0;
  double right = 1.0;
  for (;;) {
    std::uniform_real_distribution<> runif(left, right);
    double const new_t = runif(engine);
    double const new_x = std::log(1.0 / new_t - 1.0);
    double const new_y = f(new_x) / (new_t * (1 - new_t));

    if ((new_y > slice) || (t == new_t)) {
      t = new_t;
      x = new_x;
      y = new_y;
      break;
    }

    if (t < new_t) {
      right = new_t;
    } else {
      left = new_t;
    }
  }

  // record x
  std::cout << x << '\n';
}

指数型分布族

p(\bar{x}) \propto e^{f(\bar{x})}

で表される指数型分布族の確率分布からスライスサンプリングする。

一様乱数 $\xi \sim [0, 1]$ を使うと

\begin{aligned}
                     u &= e^{f(\bar{x})}       \xi \\
\Leftrightarrow \log u &=    f(\bar{x}) + \log \xi
\end{aligned}

となることから、$e^{f(\bar{y})} > u$ となる $\bar{y}$ から一様サンプリングする代わりに、$f(\bar{y}) > \log u$ となる $\bar{y}$ から一様サンプリングすればよい。

備考

任意の分布族でも $\log f(x)$ と対数をとって計算するのが普通みたい。

参考文献

  1. Slice sampling
  2. Easy unbounded slice sampling