また、一次元正規分布の累積分布関数(CDF; Cumulative Distribution Function)\( N(x)\)は、以下の形で表わされる。
\begin{eqnarray*}
N(x)&=&\int_{-\infty}^{x}f(z)dz\\
&=&\int_{-\infty}^{x}\frac {1}{\sqrt {2\pi \sigma ^{2}}}\exp\left(-{\frac {(z-\mu )^{2}}{2\sigma ^{2}}}\right)dz
\end{eqnarray*}

これを多次元に拡張する。

平均ベクトル\( \mathbf {\mu}\)、分散共分散行列\( {\Sigma}\)の多変量(\( d\)次元)正規分布の累積密度関数\( f_d(\mathbf　{x})\)は、以下の形で表される。

\begin{eqnarray*}
{\displaystyle f_d(\mathbf {x} )={\frac {1}{({\sqrt {2\pi }})^{d}{\sqrt {|\Sigma |}}}}\exp \!\left(-{\frac {1}{2}}(\mathbf {x} -{\boldsymbol {\mu }})^{\mathrm {T} }\,\Sigma ^{-1}(\mathbf {x} -{\boldsymbol {\mu }})\right)}
\end{eqnarray*}

また、多変量正規分布の累積分布関数\( N_d(\mathbf {x})\)は、以下の形で表される。
\begin{eqnarray*}
&&N_d(\mathbf {x})=\int_{-\infty}^{x_1}\cdots\int_{-\infty}^{x_d}f_d(\mathbf {z})dz_1\cdots dz_d\\
&&=\int_{-\infty}^{x}\cdots\int_{-\infty}^{x_d}{\frac {1}{({\sqrt {2\pi }})^{d}{\sqrt {|\Sigma |}}}}\exp \!\left(-{\frac {1}{2}}(\mathbf {z} -{\boldsymbol {\mu }})^{\mathrm {T} }\,\Sigma ^{-1}(\mathbf {z} -{\boldsymbol {\mu }})\right)dz_1\cdots dz_d
\end{eqnarray*}

多変量正規分布は、正規分布に従う確率変数\( \tilde{z}_1,\cdots,\tilde{z}_d\)が、それぞれ同時に\( x_1,\cdots,x_d\)以下となる確率を示しており、自然科学はもとより金融工学などの分野においても現れる重要なものである。

多変量正規分布のうち、特に相関係数\( \rho\)の二次元の標準正規分布について書けば、平均ベクトル\(\mathbf {\mu}=\left(0,0\right) \)、分散共分散行列\( \mathbf {\Sigma}=\left(\begin{array}{cc}
1 & \rho\\
\rho & 1
\end{array}\right)\)であるから、
\begin{eqnarray*}
f(x_1,x_2)&=&\frac {1}{\sqrt {2\pi(1-\rho^2)}}\exp\left(-\frac {1}{2\sqrt{1-\rho^2}}\left(x_1^2-2\rho x_1x_2+x_2^2\right)\right)
\end{eqnarray*}
\begin{eqnarray*}
&&N_2(x_1,x_2)=\int_{-\infty}^{x_1}\int_{-\infty}^{x_2}f_2(\mathbf {z})dz_1dz_2\\
&&=\int_{-\infty}^{x_1}\int_{-\infty}^{x_2}\frac {1}{\sqrt {2\pi(1-\rho^2)}}\exp\left(-\frac {1}{2\sqrt{1-\rho^2}}\left(z_1^2-2\rho z_1z_2+z_2^2\right)\right)dz_1dz_2
\end{eqnarray*}
と書ける。

SciPy.Stats.MVNの関数を用いた二次元正規分布の累積密度関数(CDF)の計算

Pythonのモジュールの一つに、科学技術計算に用いられるSciPyがある。

SciPyにはいくつかのサブモジュールが含まれ、SciPy.Statsもその一つである。

本記事で参照するのはそのさらに下位のモジュールSciPy.Stats.MVNである。MVNはMulti-Variate Normalの略と思われる。

SciPy.Stats.MVNはscipy.statsのガイドページ[1]やリファレンスマニュアル[2]に記載がないため、参考文献の[3]を参考にしている。

これによれば、多変量正規分布の平均ベクトルと分散共分散行列を指定することで多変量正規分布の累積密度関数を扱うことができ、確率計算の積分範囲を与えることで具体的な確率が計算できる。

以下ではこのscupy.stats.mvnを用いて二次元正規分布の累積密度関数を計算する。

まず、モジュールをインポートする。

from scipy.stats import mvn
import numpy as np

次に、平均ベクトルと分散共分散行列を指定する。

ここでは相関係数\( \rho=0.3\)の二次元標準正規分布を仮定する。

mu = np.array([0, 0])
S = np.array([[1,0.3],[0.3,1]])

最後に確率計算のための積分範囲を指定し、確率を出力する*1。

low=np.array([-10000,-10000]) #負の無限大の代用
upp = np.array([0.1, 0.1])
p,i=mvn.mvnun(low,upp,mu,S)

結果は次の通り。

>>> print(p)
0.33948562042391645

mvn.mvnun関数を用いれば、多変量の正規分布を用いた確率が計算できる。

確率計算の下限を負の大きな値とした上で、上限の値を変数\( \mathbf{x}\)と考えることで、累積分布関数として用いることが可能である。

近似計算（ハル『フィナンシャルエンジニアリング』の方法、Drezner1978法）

参考文献[4]の940ページに、二次元正規分布の累積分布関数の計算方法が記されたテクニカルノートのリンクが載っている[5]。

ここに示された原論文[6]によれば、二次元正規分布の累積分布関数は次のように計算する。

まず、ベースとなる関数を定義する。

相関係数\( \rho\)の二次元正規分布の累積分布ベース関数\( M_{base}(a,b;\rho)\)を
\begin{eqnarray*}
M_{base}(a,b;\rho)=\frac{\sqrt{1-\rho^2}}{\pi}\sum_{i,j=1}^{4}A_iA_jf(B_i,B_j)
\end{eqnarray*} ただし、
\begin{eqnarray*}
f(x,y)=\exp\left(a'(2x-a')+b'(2y-b')+2\rho(x-a')(y-b')\right)\\
a'=\frac{a}{\sqrt{2(1-\rho^2}},~b'=\frac{b}{\sqrt{2(1-\rho^2}}
\end{eqnarray*}
\begin{eqnarray*}
A1 = 0.3253030,~A2 = 0.4211071,~A3 = 0.1334425,~A4 = 0.006374323\\
B1 = 0.1337764,~B2 = 0.6243247,~B3 = 1.3425378,~B4 = 2.2626645
\end{eqnarray*}
で表す。

そして、引数\( a,b,\rho\)の符号によって、以下のように二次元正規分布の累積分布関数\( M(a,b;\rho)\)を定義する。

\( a,b,\rho<0\)のとき、\begin{eqnarray*}M(a,b;\rho)=M_{base}(a,b;\rho)\end{eqnarray*}
\( a<0,b>0,\rho>0\)のとき、\begin{eqnarray*}M(a,b;\rho)=N(a)-M_{base}(a,-b;-\rho)\end{eqnarray*}
\( a>0,b<0,\rho>0\)のとき、\begin{eqnarray*}M(a,b;\rho)=N(b)-M_{base}(-a,b;-\rho)\end{eqnarray*}
\( a>0,b>0,\rho<0\)のとき、\begin{eqnarray*}M(a,b;\rho)=N(a)+N(b)-1+M_{base}(-a,-b;\rho)\end{eqnarray*}
\( a\times b\times \rho>0\)のとき、\begin{eqnarray*}M(a,b;\rho)=M_{base}(a,0;\rho_1)+M_{base}(b,0;\rho_2)-\delta\end{eqnarray*} ただし
\begin{eqnarray*}\rho_1=\frac{(\rho a-b)sign(a)}{\sqrt{a^2-2\rho a b+b^2}},~\rho_2=\frac{(\rho b-a)sign(b)}{\sqrt{a^2-2\rho a b+b^2}},~\delta=\frac{1-sign(a)sign(b)}{4}\end{eqnarray*}

これをPythonで記述すれば、以下のようになる。

まず必要なモジュールをインポートする。

import numpy as np
from scipy.stats import norm
from scipy import sqrt,sign,exp,pi

そして加算記号の計算を行う。

また、引数の符号により計算結果を変える部分は、if文を用いて条件分岐させている。

def M(a,b,c):
    def M_base(a,b,c):
        A=np.array([0.3253030,0.4211071,0.1334425,0.006374323])
        B=np.array([0.1337764,0.6243247,1.3425378,2.2626645])
        def f_M(B_1,B_2):
            a_prime=a/sqrt(2*(1-c**2))
            b_prime=b/sqrt(2*(1-c**2))
            return exp(a_prime*(2*B_1-a_prime)\
            +b_prime*(2*B_2-b_prime)+\
            2*c*(B_1-a_prime)*(B_2-b_prime))
        SUM=0
        for n in range(4):
            for m in range(4):
                SUM=SUM+A[n]*A[m]*f_M(B[n],B[m])
        return sqrt(1-c**2)/pi*SUM
    if a<0 and b<0 and c<0:
        return M_base(a,b,c)
    elif a<0 and b>0 and c>0:
        return norm.cdf(x=a, loc=0, scale=1)-M_base(a,-b,-c)
    elif a>0 and b<0 and c>0:
        return norm.cdf(x=b, loc=0, scale=1)-M_base(-a,b,-c)
    elif a>0 and b>0 and c<0:
        return norm.cdf(x=a, loc=0, scale=1)\
        +norm.cdf(x=b, loc=0, scale=1)\
        -1+M_base(-a,-b,c)
    else:
        c1=(a*c-b)*sign(a)/sqrt(a**2-2*c*a*b+b**2)
        c2=(b*c-a)*sign(b)/sqrt(a**2-2*c*a*b+b**2)
        delta=(1-sign(a)*sign(b))/4
        return M_base(a,0,c1)+M_base(b,0,c2)-delta

前節の結果と同様の条件で計算してみる。

>>> M(0.1,0.1,0.3)
0.33948482307482419

SciPy.Stats.MVNの関数を用いた場合の計算と、小数第5位まで一致している。小数第6位以下の差は、計算アルゴリズムの違いであろう。

失敗例１：NumPy.Ramdom/SciPy.Statsを用いた二次元正規分布の累積密度関数(CDF)の計算

前節でも用いたNumPyおよびSciPyについて、失敗例を取り上げる。

Pythonのモジュールの一つであり、数値計算に用いられるNumyは、確率・統計に関する関数を多数備えている。

また、SciPyはNumPyの拡張版ともいうべきモジュールで、文字通り科学技術計算のための関数やサブモジュールを含む。

Numpyのrandom.gaussやrandom.multivariate_normalを使えば正規分布に従う確率変数を生成することが出来るが、確率密度関数(PDF)や累積分布関数(CDF)を指定することは出来ない。

また、SciPyのstats.norm.pdfやstats.multivariate_normalを使えば正規分布の確率密度関数(PDF)が得られるほか、stats.norm.cdfによって一次元正規分布の累積分布関数(CDF)を扱うことが出来るが、不幸にもstats.multivariate_normal.cdf、すなわち多変量正規分布の累積分布関数(CFD)は存在しない。

従って、上記方法では多変量正規分布の累積分布関数を計算することは出来ない。

失敗例２：Python/SymPyを用いた二次元正規分布の累積密度関数(CDF)の計算

SymPyとは、Pythonで数式処理を行うためのモジュールである。

SymPyを使えば（数学用語としての）変数を用いて、微積分や行列計算が行える。

この機能を利用して、まず多次元の確率密度関数を定義し、それを積分することで、多次元の累積分布関数を計算することが考えられる。

しかし、当該方法では、最終的な積分計算でエラーが起こり、累積密度関数を用いた確率の計算が出来ない。

以下に例を示すが、コーディングに当たっては上述の二次元正規分布を前提とする。

まずSymPyモジュールから関数を包括的にインポートする。

from sympy import *

次に平均ベクトルと分散共分散行列を定義する。

SymPyは行列計算を行うことが出来るため、Matrix形式で定義する。

また、分散共分散行列において、相関係数\( \rho=0.3\)としている。

m=Matrix([0,0])                             #縦ベクトルとして定義される
s_rho=Matrix([[1,0.3],[0.3,1]])

計算に必要な行列式と逆行列は、それぞれ行列オブジェクトのあとに.det()と.inv()をつけて計算する。

d=s_rho.det()        #行列式
s_inv=s_rho.inv()   #逆行列

念のため、元の行列と逆行列を書けて単位行列になっているか確かめる。

>>> s_rho*s_inv
Matrix([
[1.0,   0],
[  0, 1.0]])

二次元の確率密度関数pdfを定義する。

そのためにまず、変数\( x_1,~x_2\)をSymPy上で認識する変数として宣言する。

x1, x2 = symbols("x1 x2")
x=Matrix([x1,x2])
pdf=exp((-(x-m).T*s_inv*(x-m)/2))/(2*pi*sqrt(d))

二次元の累積分布関数を定義する。

積分計算の上限を\( z_1,z_2\)としている。

z1,z2=symbols("z1 z2")
cdf=integrate(integrate(pdf,(x1,-oo,z1)),(x2,-oo,z2))

ここでエラーとなる。

この方法でcdfを\( z_1,z_2\)の関数として定義したいが、私の環境では不可能であった。

ちなみに確率密度関数を全区間、すなわち正負の無限大を積分区間として積分すると、全確率として値は\( 1\)となる。

この計算はここまでの方法で可能であり、

>>> integrate(integrate(pdf,(x1,-oo,oo)),(x2,-oo,oo))
Matrix([[1.0]])

という結果が得られ、間違いはなさそうである。

しかし積分区間の上限を変数とした累積分布関数は、SymPyを使った方法では失敗した。

まとめ

本記事では、二次元正規分布の累積分布関数を与える2つの方法を示した。

一つはSciPyのStats.MVNモジュールを用いる方法であり、このモジュールに含まれる関数を用いて累積分布関数を作成した。

もう一つは『フィナンシャルエンジニアリング』に記載されたDrezner(1978)の方法であり、正規分布の近似式を用いる方法である。

これら方法のほかにも累積分布関数を作成できそうな方法は存在するが、その中から失敗例として2つ例を示した。