KL ダイバージェンスと最尤推定

KL ダイバージェンスと最尤推定#

パラメータの最尤推定

  • 入力

    • 真の確率分布 \(p_\ast (x)\) から生成された \(N\) 個のサンプル \(\{x^{(1)},x^{(2)},\ldots,x^{(N)}\}\)

    • パラメータ \(\theta\) で調整できる確率分布 \(p_\theta(x)\)

  • 出力

    • \(p_\theta(x)\) が \(p_\ast(x)\) に最も近づくような \(\theta\) の値

最尤推定では，以下の対数尤度を目的関数とする．

\[ \log \prod_{n = 1}^N p_\theta (x^{(n)}) = \sum_{n = 1}^N \log p_\theta(x^{(n)}) \]

この対数尤度を最大化するパラメータを以下のように表記する．

\[ \hat{\theta} = \underset{\theta}{\mathrm{argmax}} \sum_{n = 1}^N \log p_\theta(x^{(n)}) \]

ここで，「\(p_\theta(x)\) を \(p_\ast(x)\) に近づける」は， 「\(p_\theta(x)\) と \(p_\ast(x)\) の KL ダイバージェンスを最小にする」と言い換えられる．

すなわち，

\[ D_\mathrm{KL}(p_\ast || p_\theta) = \int p_\ast(x) \log \frac{p_\ast(x)}{p_\theta(x)} dx \]

を最小にするような \(\theta\) を求められればよい．

しかし，\(p_\ast\) が未知である以上，実際に積分計算を行うことはできない…

そこで，モンテカルロ法を用いて近似的に求めることを考える．

モンテカルロ法を用いた期待値の近似#

確率分布 \(p(x)\) に対する任意の関数 \(f(x)\) の期待値 \(\mathbb{E}_{p(x)}[f(x)]\) を以下のように定義する．

\[ \mathbb{E}_{p(x)}[f(x)] = \int p(x)f(x) dx \]

近似のアルゴリズム#

  1. 確率分布 \(p(x)\) からサンプル \(\{x^{(1)}, x^{(2)}, \ldots, x^{(N)}\}\) を生成する

  2. 各データ \(x^{(i)}\) における \(f(x^{(i)})\) を求め，その平均を計算する．

実験してみる#

\[\begin{split} \begin{align*} p(x) &= \frac{1}{\sqrt{2\pi}} \exp\left(-\frac{x^2}{2}\right)\\[10pt] f(x) &= x\\[10pt] g(x) &= x + 1 \end{align*} \end{split}\]

積分結果#

解析的に計算すると，

\[\begin{split} \begin{align*} \mathbb{E}_{p(x)}[f(x)] &= \int_{-\infty}^\infty p(x)f(x) ~dx\\[10pt] &= \int_{-\infty}^\infty \frac{x}{\sqrt{2\pi}} e^{-x^2 / 2} ~dx\\[10pt] &= -\frac{1}{\sqrt{2\pi}} \left[ e^{-x^2 / 2} \right]_{-\infty}^\infty\\[10pt] &= 0\\[20pt] \mathbb{E}_{p(x)}[g(x)] &= \int_{-\infty}^\infty p(x)g(x) ~dx\\[10pt] &= \int_{-\infty}^\infty \frac{x}{\sqrt{2\pi}} e^{-x^2 / 2} ~dx + \int_{-\infty}^\infty \frac{1}{\sqrt{2\pi}} e^{-x^2 / 2} ~dx\\[10pt] &= 0 + 1 = 1 \end{align*} \end{split}\]
:dep rand = "0.8.5"
:dep rand_distr = "0.4.3"

use rand::prelude::*;
use rand_distr::Normal;

// 乱数生成器を初期化
let mut rng = thread_rng();

// 標準正規分布
let std_normal = Normal::new(0.0, 1.0).unwrap();

// サンプリング回数
const NUM: usize = 50000;

fn f(x: f64) -> f64 {
    x
}

fn g(x: f64) -> f64 {
    x + 1.0
}

// fの期待値
let exp_f = std_normal
    .sample_iter(&mut rng)
    .take(NUM)
    .map(f)
    .sum::<f64>() / NUM as f64;

exp_f

-0.0005952278388999247

// gの期待値
let exp_g = std_normal
    .sample_iter(&mut rng)
    .take(NUM)
    .map(g)
    .sum::<f64>() / NUM as f64;

exp_g

0.9985557292925885

モンテカルロ法なので精度は低いが，理論値である以下の値にある程度近づいていることがわかる．

  • \(\mathbb{E}_{p(x)}[f(x)] = 0\)

  • \(\mathbb{E}_{p(x)}[g(x)] = 1\)