Note不確実性をマネージしよう。本記事の目的は、Pythonで二次元正規分布の累積分布関数をコーディングする方法を示すことである。SciPy.Stats.MVNモジュール使って解く方法と近似計算による方法を提示するとともに、累積分布関数をコーディングする際に直面する失敗例について述べる。他人の誤ちを繰り返したくない人また、以下の記事の内容を前提とする。プログラミング言語Pythonとはなにか、WindowsにPythonをインストールする方法なお、本記事の内容は下記書籍の内容を参考にしているため、合わせて参照してほしい。f(x)\)は、以下の形で表わされる。N(x)\)は、以下の形で表わされる。^{2}}}\right)dz\end{eqnarray*}これを多次元に拡張する。{x})\)は、以下の形で表される。{x})\)は、以下の形で表される。x_1,\cdots,x_d\)以下となる確率を示しており、自然科学はもとより金融工学などの分野においても現れる重要なものである。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である。Normalの略と思われる。SciPy.Stats.MVNはscipy.statsのガイドページ[1]やリファレンスマニュアル[2]に記載がないため、参考文献の[3]を参考にしている。これによれば、多変量正規分布の平均ベクトルと分散共分散行列を指定することで多変量正規分布の累積密度関数を扱うことができ、確率計算の積分範囲を与えることで具体的な確率が計算できる。以下ではこのscupy.stats.mvnを用いて二次元正規分布の累積密度関数を計算する。モジュールをインポートする。np次に、平均ベクトルと分散共分散行列を指定する。np.array([[1,0.3],[0.3,1]])最後に確率計算のための積分範囲を指定し、確率を出力する*1。0.1])p,i=mvn.mvnun(low,upp,mu,S)結果は次の通り。print(p)0.33948562042391645mvn.mvnun関数を用いれば、多変量の正規分布を用いた確率が計算できる。\mathbf{x}\)と考えることで、累積分布関数として用いることが可能である。近似計算(ハル『フィナンシャルエンジニアリング』の方法、Drezner1978法)参考文献[4]の940ページに、二次元正規分布の累積分布関数の計算方法が記されたテクニカルノートのリンクが載っている[5]。ここに示された原論文[6]によれば、二次元正規分布の累積分布関数は次のように計算する。ベースとなる関数を定義する。2.2626645\end{eqnarray*}で表す。M(a,b;\rho)\)を定義する。b+b^2}},~\delta=\frac{1-sign(a)sign(b)}{4}\end{eqnarray*}これをPythonで記述すれば、以下のようになる。必要なモジュールをインポートする。sqrt,sign,exp,piそして加算記号の計算を行う。引数の符号により計算結果を変える部分は、if文を用いて条件分岐させている。M_base(a,0,c1)+M_base(b,0,c2)-delta前節の結果と同様の条件で計算してみる。M(0.1,0.1,0.3)0.33948482307482419SciPy.Stats.MVNの関数を用いた場合の計算と、小数第5位まで一致している。小数第6位以下の差は、計算アルゴリズムの違いであろう。前節でも用いた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)は存在しない。従って、上記方法では多変量正規分布の累積分布関数を計算することは出来ない。失敗例2:Python/SymPyを用いた二次元正規分布の累積密度関数(CDF)の計算SymPyとは、Pythonで数式処理を行うためのモジュールである。SymPyを使えば(数学用語としての)変数を用いて、微積分や行列計算が行える。この機能を利用して、まず多次元の確率密度関数を定義し、それを積分することで、多次元の累積分布関数を計算することが考えられる。当該方法では、最終的な積分計算でエラーが起こり、累積密度関数を用いた確率の計算が出来ない。以下に例を示すが、コーディングに当たっては上述の二次元正規分布を前提とする。SymPyモジュールから関数を包括的にインポートする。*次に平均ベクトルと分散共分散行列を定義する。SymPyは行列計算を行うことが出来るため、Matrix形式で定義する。\rho=0.3\)としている。#縦ベクトルとして定義されるs_rho=Matrix([[1,0.3],[0.3,1]])計算に必要な行列式と逆行列は、それぞれ行列オブジェクトのあとに.det()と.inv()をつけて計算する。念のため、元の行列と逆行列を書けて単位行列になっているか確かめる。1.0]])二次元の確率密度関数pdfを定義する。x_1,~x_2\)をSymPy上で認識する変数として宣言する。x2")x=Matrix([x1,x2])pdf=exp((-(x-m).T*s_inv*(x-m)/2))/(2*pi*sqrt(d))二次元の累積分布関数を定義する。z_1,z_2\)としている。z2")cdf=integrate(integrate(pdf,(x1,-oo,z1)),(x2,-oo,z2))ここでエラーとなる。z_1,z_2\)の関数として定義したいが、私の環境では不可能であった。integrate(integrate(pdf,(x1,-oo,oo)),(x2,-oo,oo))Matrix([[1.0]])という結果が得られ、間違いはなさそうである。積分区間の上限を変数とした累積分布関数は、SymPyを使った方法では失敗した。まとめ本記事では、二次元正規分布の累積分布関数を与える2つの方法を示した。一つはSciPyのStats.MVNモジュールを用いる方法であり、このモジュールに含まれる関数を用いて累積分布関数を作成した。一つは『フィナンシャルエンジニアリング』に記載されたDrezner(1978)の方法であり、正規分布の近似式を用いる方法である。これら方法のほかにも累積分布関数を作成できそうな方法は存在するが、その中から失敗例として2つ例を示した。楽天市場で購入Yahooショッピングで購入7netで購入*1積分範囲の下限lowは負の無限大としたいところである。numpyには無限大を示すinfが定義されているためこれを使いたいが、mvnモジュールではinfを扱えないようである。仕方がないので負の大きな数をlowベクトルの要素として代用している。Next>人気記事RとRStudioをアップデートする方法。最近バージョンアップしていない人向け本記事ではRとRstudioをアップデートする方法について述べる。Rについて最近バージョンアップをしていない方(筆者は3.0.2からのアップデートであった)を対象にする。R本体をアップデートした後、RStudioをアップデートする方法を述べる。指定されているテキスト・教科書・参考書についてまとめている。その中からこのサイトおすすめのテキストを、パターン別に紹介する。本記事の目的は、Pythonで方程式を解くプログラムを作成することである。SciPyモジュールの関数を使って解く方法と、方程式の数値解法であるニュートン法、および二分法を素直にコーディングする方法の3種類を提示する。本記事では例として、以下...伊藤の公式を直感的に理解する(追記:ブラック・ショールズモデル)この記事では確率解析の重要な公式である「伊藤の公式」について述べる。伊藤の公式は確率過程の変換公式であり、日本の数学者伊藤清が創始した確率解析学の基礎となっている。この記事では直感的な理解を目指し、数学的厳密性にはこだわらない。online、Pythonオンラインを使う方法)本記事では、ブラウザ上で動作するPythonのプログラミング環境について紹介する。本来プログラミングは、自分が使うPCでプログラミング言語を動かすための「環境」を構築する必要がある。