はじめに
今回は正規分布関数(ガウス関数)の数値積分をpythonで試みた．
その際に3つの手法を比較したのでそちらを記す．
手法は「数値積分法」，「scipy」，「monte-carlo法」の3つ．
ソースコードはこちら
今回もそうであるが，実際に積分計算を行いたい式は解析的な形の積分関数を求めることができないことが一般的であるから，今回はその典型例である正規分布関数を用いた．
手法概要
今回計算したいのは以下の式：
\int_{lower}^{upper} \frac{1}{\sqrt{2\pi}\sigma}\exp\Biggl ({-\frac{(x-\mu)^2}{2\sigma^2}}\Biggl) dx
そして，引数は以下の通り：
lower(=a): 積分の下端 \\
upper(=b): 積分の上端 \\
\mu: 分布の平均値 \\
\sigma: 分布の標準偏差 \\
また，scipyを除いた他の2手法に関しては積分の評価精度に関わる変数 $n$ が存在する．
数値積分法
ここでいう数値積分法とはいわゆる区分求積法のことである．Monte-Carlo法は後述する．つまり，積分値を
\frac{b-a}{n} \sum_{k=1}^{n} f\Bigl(a + \frac{b-a}{n}k\Bigl)
によって，近似する方法である．$n$ の値が十分に大きい場合，この値は真の値に近づく．近似は以下のように行っている．$\Bigl(\Delta x=\frac{b-a}{n} \ll 1, a_k = a+k\Delta x\Bigl)$
\int_{a_k}^{a_{k + 1}}f(x)dx \simeq
\int_{a_k}^{a_{k + 1}}f(a_k)dx
= f(a_k)\Delta x
ただし，$\ 0 < \theta < 1$ として．$Taylor$展開を用いると，
f(a_k＋\theta\Delta x)\Delta x - f(a_k)\Delta x = \theta f'(a_k)\Delta x^2 + o(\Delta x^3)
この値を $\theta$ で0〜1 の積分をすると，
\int_{0}^{1}\Bigl (f(a_k＋\theta\Delta x) - f(a_k)\Bigl)\Delta x d\theta = \frac{f'(a_k)}{2!} \Delta x^2 + o(\Delta x^3)
という誤差が得られるから，全体としては，
\sum_{k=1}^{n}\frac{f'(a_k)}{2!} \Delta x^2 \simeq \Delta x \int_{a}^{b}\frac{f'(x)}{2}dx=\frac{f(b)-f(a)}{2}\Delta x
程度の誤差が生じる．つまり，定義域をどの程度細かく区切るかによって，積分の値の精度に差が出る．もちろん，プログラム上では丸め誤差も存在するため，限界ぎりぎりの桁数では精度の逆転現象が起こる可能性がある．この誤差を減らすために矩形を台形に変更して数値積分をする方法もある(今回はそのネタではないので，一旦おいておきます)．
追記
台形公式では，$a_k$ と $a_{k+1}$ の間を線形補間する方法です．そして，この手法を使用することによって，上の $\frac{f(b)-f(a)}{2}\Delta x$ の項が消えて，計算精度の桁のオーダーがワンランク上( $ \Delta x^2 $ )に上がることになる．
そのコードが下のもの．シンプルですね．
def naive_integral(lower, upper, mu, sgm, dx_div):
dx = float(1 / dx_div * (upper - lower))
x = np.linspace(lower, upper, int(dx_div))
c1 = 1.0 / np.sqrt(2 * np.pi) / sgm
idx = - 0.5 * ((x - mu) / sgm) ** 2
print("value:",np.sum(dx * c1 * np.exp(idx)))
scipy(誤差関数)
正規分布を積分した関数系は解析的に求めることができていないが，その関数を形容する単語は存在する．それが誤差関数．
誤差関数は正規分布関数の累積分布関数を求めるときに使用され，
F(x)=\int_{-\infty}^{x} \frac{1}{\sqrt{2\pi}\sigma}\exp\Biggl(-\frac{(x'-\mu)^2}{2\sigma^2}\Biggl) dx'
= \frac{1}{2}\Biggl(1+{\rm erf }\ \frac{x-\mu}{\sqrt{2}\sigma} \Biggl)
という形で累積分布関数を表したときの，${\rm erf}$ が誤差関数である．scipyにはこちらの関数値を求めるライブラリが存在し，書き方は以下の通り．
def scipy_lib(lower, upper, mu, sgm, dx_div):
idx_l = (lower - mu) / np.sqrt(2) / sgm
idx_u = (upper - mu) / np.sqrt(2) / sgm
print("value:",0.5 * ( scipy.special.erf(idx_u) - scipy.special.erf(idx_l) ) )
Monte-carlo法
最後がモンテカルロ法．一番単純な考え方で例えば，サイコロの目がそれぞれ何回出るかわからないなら，サイコロ振っちゃえみたいな手法．要は実際に乱数によるサンプリングを行い，そのサンプリングのうち条件をみたすものを数え上げ，全体に対して条件を満たす確率をとったものをその事象が起こる確率として採択を行うもの．
正規分布の積分では実際に$N(\mu,\sigma^2)$ に従う乱数を $n$ 個生成し，その乱数のうち $a〜b$ の間に入ったものがどの程度の割合存在するか測定する．
こちらの測定誤差に関しては大数の法則と中心極限定理から誤差の分散は $\frac{1}{\sqrt{n}}$ のオーダーであることがわかる．
def monte_carlo(lower, upper, mu, sgm, dx_div):
samples = np.random.normal(
loc = mu,
scale = sgm,
size = int(dx_div))
print("value:",np.mean((lower <= samples) & (samples <= upper)))
結果
正規分布表の参照はこちら．
そして，本題の計算速度．(以下の計算ではすべて $\mu=0,\sigma = 1$)
ちなみに，$\pm \sigma$ なら確率 0.68269，$\pm 5\sigma$ なら確率 1 - 5.74E-07．
一応，速度がどれくらい変わるか見るために広めの定義域で試しました．
$a = -1, b = 1, n = 1000$
正解は0.68269
function naive_integral
value: 0.6824905826186739
elapsed time[s]: 0.00016736984252929688
function scipy_lib
value: 0.6826894921370859
elapsed time[s]: 2.9802322387695312e-05
function monte_carlo
value: 0.686
elapsed time[s]: 0.0001266002655029297
$a = -1, b = 1, n = 1000000$
正解は0.68269
function naive_integral
value: 0.6826892933888814
elapsed time[s]: 0.03223371505737305
function scipy_lib
value: 0.6826894921370859
elapsed time[s]: 4.267692565917969e-05
function monte_carlo
value: 0.682288
elapsed time[s]: 0.05030488967895508
$a = -5, b = 5, n = 1000$
正解は1 - 5.74E-07
function naive_integral
value: 0.9989994420133417
elapsed time[s]: 0.0001678466796875
function scipy_lib
value: 0.9999994266968562
elapsed time[s]: 2.86102294921875e-05
function monte_carlo
value: 1.0
elapsed time[s]: 0.0001285076141357422
$a = -5, b = 5, n = 1000000$
正解は1 - 5.74E-07
function naive_integral
value: 0.9999984267122968
elapsed time[s]: 0.03359627723693848
function scipy_lib
value: 0.9999994266968562
elapsed time[s]: 5.435943603515625e-05
function monte_carlo
value: 0.999998
elapsed time[s]: 0.05029034614562988
まとめ
測定精度に関しては上で導出したとおりのオーダーでまとまっていました．こうやって，数学とプログラミングをつなげて考えるのは面白いですね．Monte-Carlo法に関しては統計物理のような分野をマスターしてからやるともっと面白そう．
速度に関してはMonte-Carlo法ではサンプル数を1000倍にすると400倍遅くなり，数値積分法では200倍遅くなった．
scipyには精度の項目が存在していないので速度は一切変わらず．ただし，精度は完璧でかつ速度も断然早いので複雑でない限りはscipyが良さそう．複雑になってきた場合に精度を強く求めない場合は少数サンプルのMonte-Carlo法で良いかもしれないが，精度をある程度必要とする場合はMonte-Carlo法ではサンプル数の増加に対して精度の安定が数値積分よりも遅いので数値積分のほうが良いかもしれない．
今回は有名なガウス関数であったため，専用のライブラリが存在したが，実際には存在しないことのほうが一般的であることを考えると，サンプル数を多めにとって精度の数値積分もしくはサンプル数少なめで良いから早く結果を求めるMonte-Carlo法，という選択を迫られることになると感じた．
追記
あとでscipyのquadという関数の存在に気づいた．
この計算速度は数値積分の項目で説明したものより断然速い．
数値積分の計算にもNumpyを利用しており，そこまでの速度差が出る理由がよくわからないが，とにかくScipyのライブラリの実行速度が速い．また，この方法であれば融通は効きづらいかもしれないが積分の関数系がわからない関数に対しても高速な積分が可能となる．
from scipy.integrate import quad
import numpy as np
def scipy_integral(f, upper, lower):
print(quad(f, upper, lower)[0])
class normal():
def __init__(self, mu, sgm):
self.mu = mu
self.sgm = sgm
def f(self):
c = 1.0 / np.sqrt(2 * np.pi) / self.sgm
def _imp(x):
idx = - 0.5 * ((x - self.mu) / self.sgm) ** 2
return c * np.exp(idx)
return _imp
mu = 0
sgm = 1
lower = -1
upper = 1
cls = normal(mu, sgm)
f = getattr(normal, "f")()
scipy_integral(f, upper, lower)