この記事の目的

現在データ分析を専門にして働いている筆者ですが、c++の勉強もかねていくつかのアルゴリズムを自分で実装して公開してみようと思います。

今回は乱数をつかって効率的に定積分を実行するモンテカルロ積分を実装します。
ゆくゆくは機械学習のアルゴリズムとかも１から実装したいですね

今回載せたコードはここにありますgithub

自己紹介

データサイエンティスト2年目です。
大学のころは物理を専攻していて、FORTRANで量子シミュレーションをしていましたが、社会人になってからはもっぱらpythonを使っています。
pythonのライブラリはとても便利ですが、私はフルスクラッチで1から実装するのが好きなのでデータサイエンス分野で触れたアルゴリズムで面白そうなものを自分で実装してみようかと思っています。

プログラミングで定積分をする方法

関数 $f(x)$を区間[a, b]で定積分するのは、言わずもがな、

\begin{align}
I = \int_a^b f(x) dx
\end{align}

ですが、これをプログラミングで計算するにはいくつかの方法があります。

区分求積法

高校で習う矩形に区切って足すあれです。

計算精度はどれだけ矩形を細かく区切るかに依存しますがその分計算コストが増加します。変数の分ループを回すのでとくに変数が多いときに計算コストが増加してしまう点がプログラウの実行時間を圧迫してしまう問題があります。

たとえば2変数で各変数の積分区間を1000この矩形に分けると、1000✕1000で1000000回被積分関数を実行する必要があります。

また、機械学習やシミュレーションでは関数を実数全体で定積分する必要があるのでその分計算コストかさんでしまいます。

ちなみに、区切る矩形の数を増やさないで計算精度をあげる工夫もいくつかあり、代表的なものには台形法があります。
1変数を愚直に実装するならこんな感じです。

quadrature.cpp

double divieded_quadrature(function<double(double)> fn, 
                           double xmin, 
                           double xmax,
                           double f){

    double s=0, x, y;
    double xrange = xmax - xmin;
    double dx = xrange / f ;
    
    for (int i = 0; i < f ; i++)
    {
        // quadrature
        x = xmin + dx * i;
        y = fn(x);
        s += y * dx; 
    }
    
    return s;
}

実際には矩形の中心(このコードだとfn(x)ではなくてfn(x + dx/2))を計算します。

モンテカルロ積分

次に乱数をつかって定積分を計算するモンテカルロ積分を紹介します
矩形に区切るところまでは区分求積法を同じです。モンテカルロ積分では計算する矩形をランダムにサンプリングします。
モンテカルロ積分のwikipedia(日本語ページなかった)

サンプリングするメリットとしては、

精度がいい感じになったところでサンプリングを打ち切れる(最初に計算量を決めなくて良い)
変数が増えた際に多重ループにしなくて良い(実装が楽)

などがあります。反面、変数が少ないときは区分求積法のほうが精度が良かったりします。

モンテカルロ積分を1変数で愚直に実装します。

quadrature.cpp


double montecarloIntegration(function<double(double)> fn, 
                             double xmin, 
                             double xmax,
                             double f){
                                
    // uniform random distribution generator
    mt19937 engin(0);
    uniform_real_distribution<> dist(xmin, xmax);
    

    double s = 0, x, y;
    for (int i = 1; i <= f; i++)
    {
        x = dist(engin);
        y = fn(x);
        s += y ; 
    }
    double dx = (xmax - xmin) / f;
    s *= dx;
    return s;
}

dxをループ回数に応じてサンプリングしたあとに決めれるのが長所ですね。
乱数の生成にはメルセンヌ・ツイスタ法を採用しています。

計算比較

ガウス関数

被積分関数としてガウス関数を使用します。今回は2変数。

\begin{align}
f(x) = \exp(- (x^2 + y^2))
\end{align}

ガウス積分は、

\begin{align}
I = \int_{- \infty}^{\infty} \int_{- \infty}^{\infty} \exp(- (x^2 + y^2))dx dy = \pi
\end{align}

で知られています

先程作成した関数を使って、

experiment.cpp

// 区分求積法
s1 = divieded_quadrature2d(gaussian2d);
// モンテカルロ積分
s2 = montecarloIntegration2d(gaussian2d);

で定積分を計算できます。実際に計算する際は積分区間を無限に取ることはできないので今回は[-100, 100]で実行します。

計算条件

積分区間[-100, 100]
区分求積法の矩形数変数ごと1000(2変数なので1000,000)
モンテカルロ積分のサンプリング数1000,000

結果

区分求積法: I = 3.14159
モンテカルロ積分 I = 3.37013

あれ、思ったより差が大きい、、
2変数くらいだと同じ区分求積法のほうが簡単に精度出せるみたいですね。
4変数（4重積分）くらいからモンテカルロ積分の恩恵を受けれるみたいです。

サンプリングに無駄が多い..?

ガウス関数はx, yが0から離れると急速に0に収束します。今回の積分では結果はほぼ[-5, 5]で決まり、それより外側の領域をいくら頑張って計算してもそこまで精度は上がりません。
したがって定積分を計算する際に値が結果に影響しやすい領域を中心にサンプリングしたほうが効率よく精度を出すことができます。

この方法を重点サンプリングといい、モンテカルロ積分を使用する際はこのメリットを享受したい場合が多いです。

重点サンプリングについては次回書こうと思います。