nobUnagaの日記

モンテカルロ法1 (モンテカルロ積分)

統計・機械学習 python

Rによるモンテカルロ法入門のメモ，その１。

イントロ

統計的推論では最適化と積分の計算，特に実用上は数値計算が必要になる。
数値積分の手法としてモンテカルロ積分がある。シンプソン法などの区分法は多次元の場合に区分を細かくしようとすると計算量が増加してしまうので適用が困難である。そのような多次元の場合に対してモンテカルロ積分は有用。

モンテカルロ 積分

モンテカルロ積分の考え方は非常に単純で，下記の期待値計算の積分を対象の確率分布からのサンプリングによる平均値で近似する方法。
例えば，f(x)が確率密度関数で，xがf(x)に従う時のh(x)の期待値を計算したいときに，

$E_f(h(x)) = \int_{x}f(x)dx$

を

$\frac{1}{m}\sum_{i=1}^{m}h(x_j)$

で近似する。

同様に考えて，定積分も対象の範囲に対する一様乱数を生成して，その時のh(x)の値を求めて平均取れば，その区間の平均値となるので，積分区間を掛けて上げれば積分値が得られる。(b-a倍するのを忘れないように)

import numpy as np
from scipy import integrate
import math

f = lambda x: math.sin(x)/x  # 被積分関数

def mc_integrate(f, a, b, n=10000):
    def scaling(x):
        '''
        np.random.randで生成する[0,1)の一様乱数を[a,b)に変換する関数
        '''
        return (b-a)*x + a

    xs = np.random.rand(n)
    xs = np.vectorize(scaling)(xs)
    return (b-a) * np.vectorize(f)(xs).sum()/n


def main():
    v = mc_integrate(f, 0, 0.0, 2.0*math.pi)
    
    # Scipyを使って検算
    v_val = integrate.quad(f, 0.0, 2.0*math.pi)[0]
    print('Scipy, v={:.5}, MC, v={:.5}'.format(v_val, v))

if __name__ == '__main__':
    main()