LoginSignup
8

More than 5 years have passed since last update.

数値積分 > QMC > Quasi Monte Calro Method

Last updated at Posted at 2015-04-12

積分の方法は色々あるが、その中でQMCという積分の紹介。

概要

10年近く前にlucilleというページで初めて知った。

Prof. Harald Niederreiterさんの「Random Number Generation and Quasi-Monte Carlo Methods (PDF, 11MB)」に詳しく掲載されている。

積分をする上でモンテカルロ法というのがある。積分するポイントを乱数で取得しながら各地点の値を足しあわせていくことで積分を行う方法である。

モンテカルロ法より少ない回数で積分精度をあげようというのがQMCという手法である。積分のオーダーについては上記のPDFを参照。

QMCではlow-discrepancy sequence(LDS)という数列を用いる。その二次元版は「Halton_sequence」という数列と認識している。

サンプルコード

以下のサンプルコードで実際に積分してみる。
Halton_sequence()で二次元のLDSを作っているが、今回は一次元のみ使っている。
なお積分範囲は[0,1]となるので、これ以外の積分の場合は[0,1]で得られるLDSを積分の希望範囲へ拡張すればいい。

sample-qmc.c
#include <stdio.h>
#include <math.h>

static void Halton_sequence(int i0, double *px0, double *py0)
{
    // configuration
    const int xbase = 2;
    const int ybase = 3;
    //
    double invxbase, invybase;
    double facx, facy;
    double x0, y0;
    int inp;

    invxbase = 1.0 / (double)xbase;
    facx = 1.0 / (double)xbase;

    invybase = 1.0 / (double)ybase;
    facy = 1.0 / (double)ybase;

    inp = i0;
    x0 = 0.0;
    while (inp > 0) {
        x0 = x0 + (inp % xbase) * invxbase;
        inp = inp / xbase;
        invxbase = invxbase * facx;
    }

    inp = i0;
    y0 = 0.0;
    while (inp > 0) {
        y0 = y0 + (inp % ybase) * invybase;
        inp = inp / ybase;
        invybase = invybase * facy;
    }

    *px0 = x0;
    *py0 = y0;
}

static double myFunc(double x_)
{
    // result [1/4]
    return (x_ * x_ * x_);
}

const static double kAnswer = 1.0/4.0;

int main(void)
{
    double xx, yy;
    int idx;
    double d0;
    int count = 0;
    double res;

    d0 = 0.0;
    for(idx=1; idx<=3000; idx++) {
        count++;
        Halton_sequence(idx, &xx, &yy);

        d0 = d0 + myFunc(xx);
        res = d0 / (double)count;
        printf("%d %f err=%.2f\%\n", idx, res, fabs(res/kAnswer - 1.0)*100);
    }
}
実行結果
1 0.125000 err=50.00%
2 0.070312 err=71.88%
3 0.187500 err=25.00%
...
300 0.245102 err=1.96%
301 0.245453 err=1.82%
302 0.244952 err=2.02%
...
1000 0.248766 err=0.49%
1001 0.248725 err=0.51%
1002 0.248517 err=0.59%
...
2998 0.249347 err=0.26%
2999 0.249531 err=0.19%
3000 0.249448 err=0.22%

ガウス積分などと比べて同じ精度では計算する点数が多くなってしまう。
QMCの良いと思う所は、点を増やすときは過去の点そのままが使えるという所。
100個で計算したけど収束していない場合、101個目から計算を続ければいい。ガウス積分やシンプソン積分だとこうはいかない。
そのため、収束精度を設定して、その精度が求まるまで計算を続けるという方法が取れそう。

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
8