分散・不偏分散・標準偏差
C++/Linuxシミュレーションプログラム
C++/Linux用のプログラムをVC++で利用する方法についてはこのウェブサイトの後半で説明している。
不偏分散
の期待値は母集団の分散
（存在するならば）に等しい（岡本安晴「データ分析のための統計学入門」おうふう、pp. 146-148 参照）
が、不偏分散の平方根（標準偏差の推定値として用いたとする）の期待値は母集団の標準偏差に等しいとはいえない。
これは、平方根を与える関数が凹関数（上に凸関数）であることからわかる（イェンセンの不等式：岡本安晴（2009）「統計学を学ぶための数学入門［下］」培風館、p. 137 参照）。
シミュレーションによって式（１）の様子を見るためのプログラム（リスト１）を作成した。
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝ リスト１ ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
#include <iostream>
#include "WHRN.h" // 乱数生成用クラス型NormalRNの宣言・定義
using namespace std;
using namespace MyLib;
int main(){
NormalRN rn; // 乱数生成用オブジェクトrn
int nsimu = 1000000; // シミュレーション（サンプリング）の回数
while(true){
int sampleSize;
cout << "Sample size = ";
cin >> sampleSize; // サンプルサイズ
if(cin.fail()) break;
if(sampleSize < 2) continue;
double *X = new double[sampleSize]; // データ設定用配列
double sum_UVar = 0.0, sum_Var = 0.0;
double sum_SqrtUVar = 0.0, sum_SqrtVar = 0.0;
for (int isimu = 0; isimu < nsimu; isimu++){
double sumX = 0.0;
for (int j = 0; j < sampleSize; j++){
X[j] = rn.normal(); // 標準正規分布
sumX += X[j];
}
double mean = sumX / sampleSize;
double sumX2 = 0.0;
for (int j = 0; j < sampleSize; j++)
sumX2 += sqr(X[j] - mean);
sum_UVar += sumX2 / (sampleSize - 1.0); // 不偏分散の和
sum_SqrtUVar += sqrt(sumX2 / (sampleSize - 1.0)); // 不偏分散の平方根の和
sum_Var += sumX2 / sampleSize; // 分散の和
sum_SqrtVar += sqrt(sumX2 / sampleSize); // 分散の平方根の和
}
cout << "mean_UVar = " << sum_UVar / nsimu << endl;
cout << "mean_Var = " << sum_Var / nsimu << endl;
cout << "mean_SqrtUVar = " << sum_SqrtUVar / nsimu << endl;
cout << "mean_SqrtVar = " << sum_SqrtVar / nsimu << endl;
delete [] X;
}
return 0;
}
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
リスト１のプログラムは、標準正規分布に従う独立なサンプルをsampleSize個採り、その不偏分散、分散およびそれらの平方根を算出することをnsimu回繰り返すものである。nsimu回の繰り返しにおける平均値をシミュレーションによる統計量の期待値の推定量としている。リスト１のプログラムの実行例を図１に示す。
図１
サンプルサイズが２（Sample size = 2）のとき、不偏分散の平均値はほぼ１で母集団の分散に等しいことが確認できる。しかし、不偏分散の平方根の平均値はmean_SqrtUVar = 0.799と母集団の標準偏差１の約８割の大きさであり、不偏推定値ではない。
リスト１のプログラムは「Sample size = 」と聞いてきたとき、２以上の整数値を入力すると、入力したサンプルサイズでのシミュレーションが行われる。「Sample size = 」に対して1以下の値を設定すると、再度入力が求められる。数値以外の文字を入力するとプログラムの実行終了となる。
リスト１でインクルードされている乱数生成用クラス型NormalRN型の宣言・定義のコードファイルWHRN.hおよびWHRN.cppは、リスト２およびリスト３に示すものである。
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝ リスト１ WHRN.h ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
#ifndef WHrn_CK
#define WHrn_CK
#include <cmath>
namespace MyLib {
using namespace std;
inline double sqr( double v ){
return v * v;
}
// Wichmann/Hill generator...the period is of the order of 10 to the 12th power.
// cf. J.E.Gentle, 2003.
// "Random Number Generation and Monte Carlo Methods, 2nd ed." Springer.
class WHrn {
public:
unsigned int x;
unsigned int y;
unsigned int z;
double u;
WHrn(): x(11111), y(22222), z(33333) {}
virtual ~WHrn(){}
WHrn(unsigned int ix, unsigned int iy, unsigned int iz):
x(ix), y(iy), z(iz) {}
double uni();
};
class NormalRN: public WHrn {
public:
NormalRN(){}
virtual ~NormalRN(){}
NormalRN(unsigned int ix, unsigned int iy, unsigned int iz){
x = ix; y = iy; z = iz;
}
// Standard Normal Distribution
double normal();
// Normal Distribution of mean m and standard deviation s
double normalMS(double m, double s);
// A Pair of Independent Random Numbers of
// the Standard Normal Distribution
void normalPair( double & n1, double & n2 );
};
} // End of namespace MyLib
#endif
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝ リスト２ WHRN.cpp ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
#include "WHRN.h"
namespace MyLib {
using namespace std;
// Wichmann/Hill generator...the period is of the order of 10 to the 12th power.
// cf. J.E.Gentle, 2003.
// "Random Number Generation and Monte Carlo Methods, 2nd ed." Springer.
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 - trunc(u);
} while (u <= 0.0);
return u;
}
// Standard Normal Distribution
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;
}
// Normal Distribution of mean m and standard deviation s
double NormalRN::normalMS(double m, double s){
double v = s * normal() + m;
return v;
}
// A Pair of Independent Random Numbers of
// the Standard Normal Distribution
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;
}
} // End of namespace MyLib
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
クラス型NormalRNはWHrn型の子孫である。WHrn型はWhichmann/Hill generator（岡本安晴（2010）「大学生のための心理学VC++プログラミング入門」勁草書房、pp. 118-123 参照）によって一様乱数を生成するものである。
リスト１～３のコードファイルは圧縮ファイルPrgSimu.zipとしてまとめた。ファイル名PrgSimu.zipをクリックしてダウンロードすることができる。
Visual
C++で利用する場合
上のC++/Linux用のプログラムをVisual C++（VC++; Visual Studio）で利用するときは、以下のようにする。要点は、
（１） 「全般｜空のプロジェクト」としてプロジェクトの新規作成を行う。
（２） プログラム中の日本語はLinux用の文字コードであるので削除する。日本語が必要であるときは、Linuxで作成された日本語を削除後、改めてVisual C++のエディタにおいて日本語で入力する。
である。
まず、プロジェクトの新規作成を以下のように行う。
Visual Studioを起動し、メニュ「ファイル｜新しいプロジェクト」を選ぶ（図２）。
図２
図３のダイアログウィンドウが表示される。
図３
左のペインから「Visual C++|全般」をクリックして表示される中央のペインの項目から「空のプロジェクト」を選ぶ。下の欄の「名前」「場所」などの設定は、他のVisual C++のプロジェクトの場合と同様である（岡本安晴「大学生のための心理学VC++プログラミング入門」勁草書房、第０章「はじめに」参照）。右下の「OK」ボタンのクリックで図４の画面になる。図３の設定の下では、ソリューション名「Project3」、プロジェクト名「Project3」であることに注意。
図４
ソリューション名「Project3」の下のプロジェクト「Project3」のフォルダを開くと図５のようになっている。
図５
このフォルダ内に、C++/Linuxのプログラムファイル（上の圧縮ファイルPrgSimu.zipを解凍したフォルダにある）をコピーする（図６）。
図６
図６の状態で、Visual Studioのメニュ「プロジェクト｜既存項目の追加」を選ぶ（図７）。
図７
図８の「既存項目の追加」ダイアログウィンドウが表示される。
図８
コピーしたファイル「SimuVarSD.cpp」、「WHRN.cpp」、「WHRN.h」を選んで、「追加」ボタンをクリックする。ソリューションエクスプローラに追加したファイルが表示されるようになる（図９）。
図９
プロジェクトのフォルダーにファイルをコピーしただけではビルドの対象にならない。上の「既存項目の追加」操作を行って、ソリューションエクスプローラに表示されるようにしておく。
ソリューションエクスプローラ内のファイル名をダブルクリックすると、そのファイルがエディで表示される（図１０）。
図１０
図１０の状態でビルドすると図１１の注意が表示される。
図１１
これは、日本語がLinuxの下で入力されたものであることによる。Visual Studioのエディタで表示はできるが、ビルドではLinuxでの日本語は問題を引き起こす。図１１における「warning～」の行の次の行の「error～」は、このLinuxでの日本語が原因である。日本語の部分を全て削除する（図１２）。
図１２
図１２の状態でビルドすると図１３のようにビルドが成功する。
図１３
F５キーを押してプログラムを実行すると図１４のように正常に実行される。
図１４
なお、標準C++準拠プログラムをVisual C++（Visual Studio）で新規に作成する方法については＜Visual C++ (Visual Studio)による標準C++準拠プログラミング＞で説明している。
注意
上のVisual C++による実行例は、Visual Studio Express 2013 Previewによるものである。Visual Studio 2012でビルドすると
「trunc識別子が見つかりません」
というようなエラーが出る。このときは、WHRN.cppファイル内の
「trunc(u)」
を
「int(u)」
に変えてビルドするとよい。