【機械学習】モンテカルロ積分を理解する

モンテカルロ積分は、モンテカルロ法を用いた積分計算を近似する手法です。

積分計算を近似して数値計算的に求める方法として、シンプソン法や台形公式などがありますが、モンテカルロ積分もこのような近似手法の1つです。

モンテカルロ法は、機械学習においては、ベイジアンモデリングにおいて、MCMC等の手法において、ある確率変数の期待値を導出する際に利用されます。

MCMC法を学ぶ際に、モンテカルロ法について理解する必要がありますが、一方で、モンテカルロ法自体もなかなか理解するのが大変なので、まずはモンテカルロ法の1つの事例であるモンテカル積分を理解していく流れで本記事を解説します。

モンテカルロ積分を理解すると、MCMCなどの手法で、確率変数の期待値をモンテカルロ法を用いて計算する流れが理解できると思います。

モンテカルロ積分の概要

モンテカルロ積分は、定積分を近似して数値計算を行うための手法です。今回は、次のような関数$f(x)$を、区間[a, b]で積分することを考えます。

ここで、モンテカルロ積分では、任意の関数$f(x)$を任意の区間で近似することができます。

そして、求めたい定積分を次の$I$で定義します。

\begin{equation}
I = \int_{a}^{b} f(x) dx
\end{equation}

ある確率分布$p(x)$が、連続確率分布である場合、関数$f(x)$は期待値の定義から次のようになります。

\begin{equation}
E[f(x)] = \int p(x) f(x) dx
\end{equation}

ここで、関数$f(x)$は、確率変数の実現値$x$における任意の関数系を指します。$f(x) = x$だと平均になりますね。この辺りがよくわからない人は、期待値の定義をご覧ください。

モンテカルロ積分の手順

モンテカルロ積分の手順は次のようになります。

Pythonでモンテカルロ積分を実装

Pythonでモンテカルロ積分を実装していきます。

実装は次のようになります。

import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as stats

np.random.seed(100)

# M: n_sampling
# s_integ: 積分区間の始まり
# e_integ: 積分区間の終わり
def montecarlo(func, s_integ, e_integ, n_sampling=1000):
    res = []
    for i in range(100):
        x = np.random.uniform(s_integ, e_integ, n_sampling)
        res.append( (e_integ -  s_integ) * func(x).sum() / n_sampling)
        
    return np.array(res).mean(), np.array(res).std()

モンテカルロ積分を試す

ではこちらの実装でモンテカルロ法の結果を試してみましょう。

f(x) = 2 x^2 + x

関数$f(x)$を区間[0, 1]で定積分した値をモンテカルロ積分で求めていきましょう。この関数$f(x)$の区間[0,1]における定積分の値は、

\begin{split}
\int_{0}^{1} f(x) dx &= 
\biggl[
F(x)
\biggr]_{0}^{1}  \\
&= 
\biggl[
\frac{2}{3} x^3 + \frac{1}{2}x^2
\biggr]_{0}^{1} \\
&= \frac{7}{6} = 1.1666...
\end{split}

となります。これをモンテカルロ積分を用いて求めていきましょう。

def func1(x):
    return 2 * (x ** 2) + x

m, s = montecarlo(func1, 0, 1, n_sampling=1000000)
print("estimation: {} ±  {}".format(m, s))

estimation: 1.1665560613356636 ± 0.0008882157667182995

となり、ほとんど正確に求めることができました。

続いて正規分布の場合でもやってみます。

# 正規分布の場合
def norm_func(x):
    return stats.norm(0, 1).pdf(x)

gt = stats.norm(0, 1).cdf(2) - stats.norm.cdf(-2)
print(f"ground truth: ", gt)

m, s = montecarlo(norm_func, -2, 2, n_sampling=100000)
print("estimation: {} ±  {}".format(m, s))

ground truth: 0.9544997361036416
estimation: 0.9545133102756458 ± 0.0014196682783825184

となり、かなり正確に定積分を求めることができました。