正規分布の確率計算
C++/Linuxプログラミング
標準正規分布
の累積確率分布
を求めるプログラムと、
が与えられたときに（２）式において対応する
を求めるプログラムをC++/Linuxで作成した。
式（２）において が与えられたときに対応する を求めるプログラム
このプログラムは、式（２）における積分をガウス・ルジャンドルの積分公式を用いた適応的方法（岡本安晴（2009）「統計学を学ぶための数学入門［下］」培風館、pp.
117-118）によって行うものである。ソースリストをリスト１に示す。
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝ リスト１ CalcNormalCum.cpp ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
#include <iostream>
#include <iomanip>
#include <cmath>
#include "CalcIntgrt.h" // 積分用クラス型 AdaptiveGL の宣言
using namespace std;
using namespace myLib;
double knormal(
double z ){
return
exp(-0.5 * z * z);
}
double Cum_Normal(
double z ) {
double
p;
double
pi = atan(1.0) * 4.0;
if
(z == 0.0){
p
= 0.5;
}else
if(z > 0.0){
AdaptiveGL agl; // 積分用クラス型 AdaptiveGL のオブジェクト agl
p
= 0.5 + agl.integral( 0.0, z, knormal
) / sqrt(2.0 * pi);
}
else {
AdaptiveGL agl;
p
= 0.5 - agl.integral( z, 0.0, knormal
) / sqrt(2.0 * pi);
}
return
p;
}
int main(){
while(true){
cout << "z = ";
double
z;
cin >> z;
if(cin.fail()) break;
double
p = Cum_Normal( z ); // p
= Prob(X < z), X ~ standard normal distribution
cout << "p = " << fixed << setprecision(9) << p << endl;
}
return
0;
}
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
積分をクラス型AdaptiveGLのオブジェクトaglによって行っている。このクラス型AdaptiveGLはリスト２、リスト３に示すヘッダーファイルとソースコードにより宣言・定義されている。
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝ リスト２ CalcIntgrt.h ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
#ifndef
CALCINTGRL_CK
#define CALCINTGRL_CK
namespace myLib
{
//
// ガウス・ルジャンドルの積分公式による適応的方法
// 岡本安晴（2009）「統計学を学ぶための数学入門[下]」培風館、pp. 117-118.
//
class
AdaptiveGL {
public:
double
(*f)(double);
AdaptiveGL(){}
virtual
~AdaptiveGL(){}
double
integral( double a, double b, double (*func)(double));
void calcIntegral(
double & s0, double a, double b );
double
Gauss_Legendre( double a, double b );
};
}
#endif
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝ リスト３ CalcIntgrt.cpp ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
#include "CalcIntgrt.h"
#include <cmath>
using namespace std;
using namespace myLib;
namespace myLib
{
//
// ガウス・ルジャンドルの積分公式による適応的方法
// 岡本安晴（2009）「統計学を学ぶための数学入門[下]」培風館、pp. 117-118.
//
double
AdaptiveGL::integral(double a, double b, double (*func)(double)){
f
= func;
double
v = 0.0;
calcIntegral( v, a, b );
return
v;
}
void
AdaptiveGL::calcIntegral(double
& s0, double a, double b){
double
s1 = Gauss_Legendre( a, 0.5*(a + b) );
double
s2 = Gauss_Legendre( 0.5*(a + b), b );
if
(((abs( s0 - s1 - s2 ) < (1.0e-14)*abs(s1 + s2))
||
((abs(s0)
+ abs(s1 + s2)) < 1.0e-14))
&&
(abs(b
- a) < 1.0)){
s0
= s1 + s2;
return;
}
else{
calcIntegral( s1, a, 0.5*(a + b) );
calcIntegral( s2, 0.5*(a + b), b );
s0
= s1 + s2;
return;
}
}
double
AdaptiveGL::Gauss_Legendre(double
a, double b){
double
x0[9] = {
-9.68160239507626090E-0001,
-8.36031107326635794E-0001,
-6.13371432700590397E-0001,
-3.24253423403808929E-0001,
0.00000000000000000E+0000,
3.24253423403808929E-0001,
6.13371432700590397E-0001,
8.36031107326635794E-0001,
9.68160239507626090E-0001
};
double
w[9] = {
8.12743883615744120E-0002,
1.80648160694857404E-0001,
2.60610696402935462E-0001,
3.12347077040002840E-0001,
3.30239355001259763E-0001,
3.12347077040002840E-0001,
2.60610696402935462E-0001,
1.80648160694857404E-0001,
8.12743883615744120E-0002
};
double
v1 = 0.0;
double
i_w = b - a;
double
xi;
for
(int i = 0; i < 9; i++){
xi
= (0.5 * i_w * x0[i]) +
(0.5 * (b + a));
v1
+= w[i] * f(xi);
}
v1
= 0.5 * i_w * v1;
return
v1;
}
}
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
プログラムを実行すると標準入力ストリームcinから
cin >> z;
とzの値が読み込まれるが、数値として不適切な入力（英字など）があるとcin.fail()はtrueを返すので、次のif文
if(cin.fail()) break;
によって、while文を抜け出し、プログラムの終了となる。
プログラムの実行例を図１に示す。
図１
zの値として英字qが入力されたとき、プログラムが終了していることがわかる。
プログラムのソースファイルを圧縮ファイルPrgNormalCum.zipとしてまとめた。ファイル名PrgNormalCum.zipをクリックしてダウンロードすることができる。
式（２）において が与えられたときに対応する を求めるプログラム
次式
を満たす に対する をニュートン法（岡本安晴（2009）「統計学を学ぶための数学入門［下］」培風館、pp.
82-83）によって求める。ニュートン法が収束しないときは、Bisection法に切り替える。プログラムをリスト４に示す。
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝ リスト４ calcz.cpp ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
#include <iostream>
#include <iomanip>
#include <cmath>
#include "calcNormal.h" // ニュートン法の関数
Newton の宣言
using namespace std;
using namespace calcnormal;
double dcalc_p_normal(
double z ){
double
pi = atan(1.0) * 4.0;
double
v = exp(-0.5 * z * z) / sqrt(2.0
* pi);
return
v;
}
double calc_p_normal(
double z ){
double
p = Cum_Normal( z );
return
p;
}
double calc_z_normal(
double p ){
if
((p <= 0.0) || (p >= 1.0)) {
throw
"Invalid value for p.";
}
double
z = 0.0;
Newton(
z, p, calc_p_normal, dcalc_p_normal,
1.0e-12 ); // ニュートン法による根の計算
return
z; // p
= Prob(X < z), X ~ standard normal distribution
}
int main(){
while(true){
cout << "p = ";
double
p;
cin >> p;
if(cin.fail()) break;
try{
// p
= Prob(X < z), X ~ standard normal distribution
double
z = calc_z_normal( p );
cout << "z = " << setprecision(9)
<< z << endl;
}
catch(const char * ce){
cout << ce << endl;
continue;
}
}
return
0;
}
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
p の入力として英字など不適当な入力があったときは、cin.fail()がtrueを返すので、if文
if(cin.fail()) break;
により、while文の外に出て、プログラムの実行終了となる。
p に設定された値が確率として不適切な値のときは、関数 calc_z_normal
において例外がthrowされ、main関数におけるcatchでの例外処理によってcontinue文が呼び出されて、while文の実行が次のループに移る。
ニュートン法はリスト５およびリスト６に示すコードによって行われている。
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝ リスト５ calcNormal.h ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
#ifndef
CALCNORMAL_CK
#define CALCNORMAL_CK
namespace calcnormal
{
double
Cum_Normal( double z );
//
// Newton法による根の計算
//
void
Newton( double
& x, // 根の初期値、戻り値は根 f(x)=c
double
c, // 関数値
double
(*f)( double ), // 関数
double
(*df)( double ), // 導関数
double
acc );
}
#endif
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝ リスト６ calcNormal.cpp ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
#include <cmath>
#include "CalcIntgrt.h" // ガウス・ルジャンドル積分公式用クラス型 AdaptiveGL の宣言
#include "calcNormal.h"
namespace calcnormal
{
using
namespace std;
using
namespace myLib;
using
namespace calcnormal;
double
knormal( double z ){
return
exp(-0.5 * z * z);
}
double
Cum_Normal( double z ) {
double
p;
double
pi = atan(1.0) * 4.0;
if
(z == 0.0){
p
= 0.5;
}else
if(z > 0.0){
AdaptiveGL agl; // 積分用クラス型 AdaptiveGL のオブジェクト agl
p
= 0.5 + agl.integral( 0.0, z, knormal
) / sqrt(2.0 * pi);
}
else {
AdaptiveGL agl;
p
= 0.5 - agl.integral( z, 0.0, knormal
) / sqrt(2.0 * pi);
}
return
p; // p
= Prob(X < z), X ~ standard normal distribution
}
//
// Bisection法による根の計算
//
void
Bisection( double (*f)(double),
double s,
// fが増加関数のとき
s > 0.0
double
c,
// 根に対する関数値
double
L_b, double U_b,
// 根の存在区間の下限と上限
double
& Root,
// 根 f(Root) = c
double
acc ){
if
( (f(L_b) - c) * (f(U_b) -
c) > 0.0 ){
throw "上限と下限の組み合わせが不適切です。";
}
do
{
double
m = 0.5 * (L_b + U_b);
double
v = f(m);
if
( s * (v - c) > 0.0 ) {
U_b = m;
}
else {
L_b = m;
}
}
while ( (abs(U_b - L_b)
>= acc * abs(U_b))
&&
( (abs(U_b) +
abs(L_b)) >= 1.0e-14) );
Root
= 0.5 * (L_b + U_b);
}
//
// Newton法による根の計算
//
void
Newton( double
& x, // 根の初期値、戻り値は根 f(x)=c
double
c, // 関数値
double
(*f)( double ), // 関数
double
(*df)( double ), // 導関数
double
acc ){
int NStep = 0;
int ck = 0;
double x0;
do {
NStep++;
if (NStep >
20) { // 20回以上繰り返しのとき
ck = 1;
break;
};
x0 = x;
x = x0 + (c - f(x0)) / df(x0);
} while ( (abs(x - x0) >= (acc * (abs(x) + abs(x0))))
&&
(
(abs(x) + abs(x0)) >= 1.0e-14 )
);
if (ck == 0)
return;
// Bisection法で計算
double L_b = x;
do {
L_b -= 1.0;
} while ( f(L_b)
>= c );
double U_b = x;
do {
U_b += 1.0;
} while ( f(U_b)
<= c );
Bisection( f, 1.0, c, L_b,
U_b, x, acc );
return;
}
}
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
積分のためのクラス型を宣言・定義しているソースコードファイルは、先にリスト２およびリスト３に示したものである。
プログラムの実行例を図２に示す。
図２
不適切な範囲のpの値「p = 9」が入力されたときは例外の処理が行われていること、pに英字qが入力されたときはプログラムの実行が終了していることがわかる。
プログラムファイルは圧縮ファイルPrgNormalDev.zipにまとめた。ファイル名PrgNormalDev.zipをクリックしてダウンロードすることができる。