必要になったのでこちらの記事を投稿します．

問題設定

条件と解きたい問題は以下の通り．

ある変数\ x_1,x_2,x_3\ が存在し，それぞれの変数は互いに独立でかつ\\
x_1 〜 N(\mu_1,\sigma_1^2) \\
x_2 〜 N(\mu_2,\sigma_2^2) \\
x_3 〜 N(\mu_3,\sigma_3^2) \\
という確率分布に従う． 

そのとき，x_1 < x_2 < x_3\ となる確率\ p\ を求めよ．

また，正規分布関数を $f(x|\mu,\sigma)$ ，正規分布関数の累積分布関数を$F(X \leq x|\mu,\sigma)=F(x|\mu,\sigma)~\int_{-\infty}^{x}f(x|\mu,\sigma)dx$ と定義する．

解法

2つの独立な正規分布に従う変数の大小関係は簡単に求めることができる．なぜならば，その2つの変数の差もまた正規分布に従うからである．一方で3つの独立な変数の場合はこのような求め方はできないので数値積分もしくはモンテカルロ法のどちらかをする必要がある．

$x_2$ がある値 $t$ を取るとき，$x_1$ が $t$ 未満，$x_2$ が $t$ より大きな値を取れば良い．

したがって，

p=
\int_{-\infty}^{\infty}\int_{x_2}^{\infty}\int_{-\infty}^{x_2}
f(x_1|\mu_1,\sigma_1)f(x_2|\mu_2,\sigma_2)f(x_3|\mu_3,\sigma_3)dx_1dx_2dx_3 \\
=\int_{-\infty}^{\infty}f(x_2|\mu_2,\sigma_2)F(x_2|\mu_1,\sigma_1) \bigl( 1-F(x_2|\mu_3,\sigma_3) \bigl)dx_2

のような積分をすればよい．これらの積分をモンテカルロ法もしくは数値積分によって行う．

数値積分の方法と精度のオーダーはこちらに記したとおり．

今回の計算に関しても圧倒的に数値積分の精度が高い．ただし，少ないサンプル数の際にはMonte-Carlo法が圧倒的に高速で計算可能であるため，特に1%精度を求めていない場合はMonte-Carlo法で十分であり，それ以上の精度を必要とする場合には数値積分を行う必要がある．

import time
import numpy as np
import scipy.special

def measure_time(func, mus, sgms, dx_div):
	start = time.time()
	func(mus, sgms, dx_div)
	finish = time.time()
	print("elapsed time[s]:",finish - start)

def normal_cdf(x, mu, sgm):
	idx = (x - mu) / np.sqrt(2) / sgm
	return 0.5 * (1 + scipy.special.erf(idx))

def normal_pdf(x, mu, sgm):
	c = 1.0 / np.sqrt(2 * np.pi) / sgm
	idx = - 0.5 * ((x - mu) / sgm) ** 2
	return c * np.exp(idx)

def integral(mus, sgms, dx_div):
	lower = min([mu - 3 * sgm for mu, sgm in zip(mus, sgms)])
	upper = min([mu + 3 * sgm for mu, sgm in zip(mus, sgms)])
	dx = float(1 / dx_div * (upper - lower))
	x = np.linspace(lower, upper, int(dx_div))

	F1 = normal_cdf(x, mus[0], sgms[0])
	f2 = normal_pdf(x, mus[1], sgms[1]) 	
	F3 = normal_cdf(x, mus[2], sgms[2])

	print("value:", np.sum(dx * F1 * f2 * (1 - F3)))


def monte_carlo(mus, sgms, dx_div):
	samples = [np.random.normal(loc = mu,scale = sgm,size = int(dx_div)) for mu, sgm in zip(mus, sgms)]
	print("value:",np.mean((samples[0] <= samples[1]) & (samples[1] <= samples[2])))


mus = np.array([1.1,1.2,1.3])
sgms = np.array([0.1,0.1,0.1]) 
dx_div = 1e+6

functions = [
			"integral", 
			"monte_carlo"]

for func in functions:
	print("function",func)
	measure_time(
		eval(func),
		mus, 
		sgms, 
		dx_div
		)
	print("")