２項分布の信頼区間
C++/Linuxプログラミング
２項分布
のパラメータ
について、信頼係数
の信頼区間を求める。
N回の試行においてk回当該の事象が生起したときの信頼区間
は、次式
で与えられる（岡本、1998、pp. 176-179）。
上記の方法で与えられる信頼区間と確率
との関係は、「芝・南風原（1990）、pp. 90-92」で論じられている。
上記の信頼区間を算出するC++/Linuxプログラム例をリスト１に示す。
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝ リスト１ CalcConfInvl.cpp ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
#include <iostream>
#include <iomanip>
#include <cmath>
#include "common.h"
using namespace std;
using namespace cmmn;
int main(){
cout << "信頼性係数 p(%) = ";
double p;
cin >> p;
if ((p < 50.0) || (100.0 <= p)){
cout << "信頼性係数の値が不適切です。 50% 以上、100% 未満の値を設定して下さい。" << endl;
return 0;
}
double pa = (100.0 - p)/100.0; // 有意水準
cout << "N = ";
cin >> N;
if (N < 1){
cout << "N の値が不適切です。 N は 1 以上の値を設定します。" << endl;
return 0;
}
cout << "k = ";
cin >> k;
if ((k < 0) || (k > N)){
cout << "k の値が不適切です。 0 <= k <= N であるようにします。" << endl;
return 0;
}
LogCombi = new double[N+1]; // ２項係数の対数値用配列
for (int i = 0; i <= N; i++){
if ((i == 0) || (i == N)) LogCombi[i] = 0.0;
else {
double v = 0.0;
for (int j = i; j >= 1; j--)
v += log(N + 1.0 - j) - log( j );
LogCombi[i] = v;
}
}
double LB, UB;
if (UCum_prob(0.0) > (pa/2.0)){
LB = 0.0;
Bisection( Cum_prob, -1.0, pa/2.0, 0.0, 1.0, UB, 1.0e-9 );
}else{
if (Cum_prob(1.0) > (pa/2.0)){
UB = 1.0;
Bisection( UCum_prob, 1.0, pa/2.0, 0.0, 1.0, LB, 1.0e-9 );
}else{
Bisection( Cum_prob, -1.0, pa/2.0, 0.0, 1.0, UB, 1.0e-9 );
Bisection( UCum_prob, 1.0, pa/2.0, 0.0, 1.0, LB, 1.0e-9 );
}
}
cout << p << "% 信頼区間 = [" << setprecision(5) << LB << ", " << UB << "]" << endl;
delete [] LogCombi;
return 0;
}
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
上のプログラムCalcConfInvl.cppの実行例を図１に示す。
図１
図１では、N＝２０試行の２項分布においてｋ回の頻度に対する確率ｐの９５％信頼区間が ［0.23058, 0.68472］ であることが示されている。
リスト１におけるヘッダーファイル common.h とそれに対するソースコードファイルcommon.cppはリスト１およびリスト２に示すものである。
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝ リスト１ common.h ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
#ifndef COMMON_CK
#define COMMON_Ck
namespace cmmn {
extern double p, pa;
extern int N, k;
extern double * LogCombi;
// 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 );
//
// ２項分布
//
double prob_term( double p, int N, int i );
//
// Prob(X <= k)
//
double Cum_prob( double p );
//
// Prob(X >= k)
//
double UCum_prob( double p );
}
#endif
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝ リスト２ common.cpp ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
#include <cmath>
#include "common.h"
using namespace std;
namespace cmmn {
double p, pa;
int N, k;
double * LogCombi;
// 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);
}
//
// ２項分布
//
double prob_term( double p, int N, int i ){
if (p <= 0.0){
if (i >= 1) return 0.0;
else return 1.0;
}else if (p >= 1.0){
if (i == N) return 1.0;
else return 0.0;
}else{
return exp(LogCombi[i] + (i * log(p) + (N - i) * log(1.0 - p)));
}
}
//
// Prob(X <= k)
//
double Cum_prob( double p ){
double v = 0.0;
for (int i = 0; i <= k; i++)
v += prob_term( p, N, i );
return v;
}
//
// Prob(X >= k)
//
double UCum_prob( double p ){
double v = 0.0;
for (int i = N; i >= k; i--)
v += prob_term( p, N, i );
return v;
}
}
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
上のプログラムのファイルは圧縮ファイルConfIntvlBinmlPrg.zipとしてまとめた。ファイル名ConfIntvlBinmlPrg.zipをクリックしてダウンロードすることができる。
参考文献
岡本安晴（1998）「Delphiで学ぶデータ分析法」CQ出版
芝 祐順・南風原朝和（1990）「行動科学における統計解析法」東京大学出版会