「第８章 積分と乱数」Win32コンソールアプリケーション版
＜岡本安晴「大学生のための心理学VC++プログラミング入門」＞の第８章のサンプルプログラムのWin32版をVisual C++2008で作成した。プログラムファイルは圧縮ファイルとしてまとめた。圧縮ファイルはクリックしてダウンロードしたものをマウスの右ボタンでクリックして表示されるメニュから「解凍」あるいは「展開」を選んで解凍することができる。解凍したプログラムファイルはVisual C++2010でも開くことができる。
リスト8.1-1のWin32版
圧縮ファイルsample8_1_1win32.zip
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
#include "stdafx.h"
#include <iostream>
#include <iomanip> // setprecison
の使用のため
#include <cmath> // 数学関数の使用のため
using namespace std;
//
// ガウス・ルジャンドルの積分公式による適応的方法
//
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 );
};
//
// 関数funcのaからbまでの積分を求める。
//
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;
}
//
// 標準正規分布の確率密度関数から定数係数を除いたもの
//
double normalDistriK( double
z ){
return exp(-0.5 * z * z);
}
//
// 標準正規累積確率関数:
//
double CumNormal( double
z ){
AdaptiveGL
agl;
const double PI = 4.0
* atan(1.0);
double v = 0.5;
if (z > 0.0){
v
= 0.5 + agl.integral( 0.0, z, normalDistriK ) / sqrt(2 * PI);
}else if (z < 0.0){
v
= 0.5 - agl.integral(z, 0.0, normalDistriK) / sqrt(2 * PI);
}
return v;
}
int _tmain(int argc, _TCHAR* argv[])
{
cout
<< "z = ";
double z;
cin
>> z;
//
// p
= Prob( X <= z ); X 〜標準正規分布
//
double p = CumNormal( z );
cout
<< "p = " <<
setprecision(12) << p << endl;
cout
<< endl << "何か文字を入力して終了。" << endl;
char ck;
cin
>> ck;
return 0;
}
リスト8.2-1のWin32版
圧縮ファイルsample8_2_1win32.zip
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
#include "stdafx.h"
#include <iostream>
#include <iomanip> // setw
の使用のため
#include <cmath> // 数学関数使用のため
using namespace std;
//
// Wichmann/Hill generator...周期= 約.0e12
//
class WHrn {
public:
unsigned __int32 x;
unsigned __int32 y;
unsigned __int32 z;
double u;
WHrn():
x(11111), y(22222), z(33333) {}
virtual ~WHrn(){}
WHrn(unsigned __int32 ix, unsigned __int32 iy, unsigned __int32 iz):
x(ix),
y(iy), z(iz) {}
double uni(); //
0 < uni() < 1
};
//
// 正規乱数
//
class NormalRN: public WHrn {
public:
NormalRN(){}
virtual ~NormalRN(){}
NormalRN(unsigned __int32 ix, unsigned __int32 iy, unsigned __int32 iz):
WHrn(ix, iy, iz) {}
double sqr( double x
){
return x * x;
}
//
// 標準正規乱数
//
double NormalRN::normal();
//
// 平均m、 標準偏差s の正規乱数
//
double NormalRN::normalMS(double
m, double s);
//
// 独立な標準正規乱数n1 とn2 の対
//
void
NormalRN::normalPair( double & n1, double & n2 );
};
double WHrn::uni(){
do {
x
= (171 * x) % 30269;
y
= (172 * y) % 30307;
z
= (170 * z) % 30323;
u
= (x/30269.0) + (y/30307.0) + (z/30323.0);
u
= u - int(u);
}
while (u <= 0.0);
return u;
}
//
// A
Rejection Polar Method for Normal Variates
//
double NormalRN::normal(){
double v1, v2, w;
do {
v1
= 2.0 * uni() - 1.0;
v2
= 2.0 * uni() - 1.0;
w = sqr(v1) + sqr(v2);
}
while (((w >= 1.0) || (w <= 0.0)));
double c = sqrt(-2.0 * log(w)/w);
return c * v1;
}
double NormalRN::normalMS(double m, double s){
double v = s * normal() + m;
return v;
}
//
// A
Rejection Polar Method for Normal Variates
//
void NormalRN::normalPair( double & n1, double
& n2 ){
double v1, v2, w, c;
do {
v1 = 2.0 * uni() - 1.0;
v2 = 2.0 * uni() - 1.0;
w
= sqr(v1) + sqr(v2);
}
while (((w >= 1.0) || (w <= 0.0)));
c
= sqrt(-2.0 * log(w) / w);
n1
= c * v1;
n2
= c * v2;
}
int _tmain(int argc, _TCHAR* argv[])
{
WHrn
rn;
NormalRN
rn1, rn2(1, 2, 3);
cout
<< "rn: x = " << rn.x
<< " \ty = " << rn.y
<< " \tz = " << rn.z
<< endl;
cout
<< "rn1:x = " << rn1.x
<< " \ty = " << rn1.y
<< " \tz = " << rn1.z
<< endl;
cout
<< "rn2:x = " << rn2.x
<< " \ty = " << rn2.y
<< " \tz = " << rn2.z
<< endl;
for (int i = 0; i
< 20; i++){
cout
<< setw(7) << rn.uni() << " \t" << rn1.normal();
double n1, n2;
rn2.normalPair(n1,
n2);
cout
<< " \t" << n1 << "
\t" << n2 << endl;
}
cout
<< endl << "何か文字を入力して終了。" << endl;
char ck;
cin
>> ck;
return 0;
}