２項分布の信頼区間
openSUSE13.2でのg++のインストール、日本語入力の設定の説明もこのウェブサイトの後半で行っている。
２項分布
のパラメータ
について、信頼係数
の信頼区間を求める。
N回の試行においてk回当該の事象が生起したときの信頼区間
は、次式
で与えられる（岡本、1998、pp. 176-179）。
上記の方法で与えられる信頼区間と確率
との関係は、「芝・南風原（1990）、pp. 90-92」で論じられている。
確率の計算（２）および（３）における計算の精度を上げるため、double 型ではなく、long double 型を用いる。long double 型の精度を確かめるプログラムを、このウェブサイトの後半に用意している。
式（２）、（３）によって信頼区間を求めるプログラムを以下に示す。
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝ ファイル main.cpp ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
#include
<iostream>
#include
<cmath>
#include
"common.h"
using
namespace std;
using
namespace cmmn;
int
main(){
while (true) {
L1: cout
<< "p(%) = ";
cin >> p;
if
(p <= 0.0) break;
if
((p < 50.0) || (100.0 <= p)){
goto L1;
}
L2: pa = (100.0 - p)/100.0;
cout << "N = ";
cin >> N;
if
(N < 1){
goto L2;
}
L3: cout
<< "k = ";
cin >> k;
if
((k < 0) || (k > N)){
goto L3;
}
LogCombi = new long double[N+1];
for
(int i = 0; i <= N; i++){
if
((i == 0) || (i == N)) LogCombi[i] = 0.0;
else
{
long
double v = 0.0;
for
(int j = i; j >= 1; j--)
v
+= log(N + 1.0 - j) - log( j );
LogCombi[i] = v;
}
}
long
double LB, UB;
try{
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 << "% Interval[" <<
LB << ", " << UB << "]" << endl;
}
catch
(char const * e){
cout << "Error-> " << e << endl;
}
delete
[] LogCombi;
}
return 0;
}
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝ ファイル main.cpp ここまで ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
ヘッダーファイルcommon.hとそのソースファイルcommon.cppは、以下に示す。
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝ ファイル common.h ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
#ifndef
COMMON_CK
#define COMMON_Ck
namespace cmmn
{
extern long double p, pa;
extern int
N, k;
extern long double * LogCombi;
// Bisection
void Bisection( long double (*f)(long
double),
long double s, // f is increasing -> s > 0.0
// f is decreasing -> s < 0.0
long double c, // f(Root) = c
long double L_b,
// L_b < Root
< U_b
long double U_b,
long double & Root, // Set Root so that f(Root) = c
long double acc
); // accuracy
long double prob_term(
long double p, int N, int i );
long double Cum_prob(
long double p );
long double UCum_prob(
long double p );
}
#endif
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝ common.hここまで ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝ ファイルcommon.cpp ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
#include
<sstream>
#include
<cmath>
#include
"common.h"
using
namespace std;
namespace
cmmn {
long double p,
pa;
int N, k;
long double * LogCombi;
// Bisection
void Bisection(
long double (*f)(long double),
long double s, // f is increasing -> s > 0.0
// f is decreasing -> s < 0.0
long double c, // f(Root) = c
long double L_b,
// L_b < Root
< U_b
long double U_b,
long double & Root, // Set Root so that f(Root) = c
long double acc
){ // accuracy
if
( (f(L_b) - c) * (f(U_b) -
c) > 0.0 ){
ostringstream sout;
sout << "Invalid values...L_b
= " << L_b << " U_b = "
<< U_b;
string
strmsg = sout.str();
char
* cstr = new char[strmsg.length()
+ 1];
for
(int i = 0; i < strmsg.length(); i++)
cstr[i] = strmsg[i];
cstr[strmsg.length()] = 0;
throw
cstr;
}
do
{
long
double m = 0.5 * (L_b + U_b);
long
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);
}
long double prob_term( long 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)));
}
}
long double Cum_prob( long double p ){
long
double v = 0.0;
for
(int i = 0; i <= k; i++)
v
+= prob_term( p, N, i );
return
v;
}
long double UCum_prob( long double p ){
long
double v = 0.0;
for
(int i = N; i >= k; i--)
v
+= prob_term( p, N, i );
return
v;
}
}
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝ common.cppここまで ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
上の３つのファイルを、圧縮ファイルPrgBinomial.zipにまとめた。圧縮ファイル名PrgBinomial.zipをクリックしてダウンロードすることができる。
2項分布（１）における2項係数は、そのまま計算するとオーバーフローを起こしやすい。２項係数の計算におけるオーバーフローと確率計算におけるアンダーフローを避けるため、式（１）を対数変換を利用した式
を用いて計算している。
上のプログラムをopenSUSE 13.2 で実行した例を図１に示す。
図１
コマンド g++ は、openSUSE13.2 (DVD)をインストールしただけでは使えないので、g++ のインストールを行う必要がある。g++ のインストールは、以下の手順で出来る。
まず、「YaST2モジュ…」を選ぶ（図２）。
図２
次に、「ソフトウェア管理」を選ぶ（図３）。
図３
表示された「YaST2」において、検索ページの検索欄に「gcc-c++」を設定して検索すると、右のペインに「gcc-c++」と「gcc-c++-32bit」が表示されるので、チェックを入れて選択する（図４）。
図４
上の設定の下で「了解」ボタンをクリックするとインストールが始まる。
また、openSUSE13.2のインストール後は、画面上のバーの「_ja」をクリックすると図５のように表示されるが、この状態では日本語入力はできない。
図５
日本語入力は、以下の手順でできるようになる。
まず、「設定」を選ぶ（図６）。
図６
次に、「地域と言語」を選ぶ（図７）。
図７
「入力ソース」で「＋」をクリックする（図８）。
図８
表示された「入力ソースの追加」において、「日本語」を選び（クリックし）（図９）、続いて「日本語（Mozc – JP layout）」を選んで「追加」をクリックする（図１０）。
図９
図１０
上の操作の後は、「日本語（Mozc – JP layout）」が選べる（図１１）。
図１１
「日本語（Mozc – JP layout）」を選んだ状態では、「半角／全角 漢字」キーの押下で、日本語入力（全角）、英数字（半角）の切り替えができる。
long
double 型の精度
long double 型の精度を確かめるプログラムを以下のように用意した。
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
#include
<iostream>
#include
<iomanip>
using namespace
std;
int
main(){
long
double v = 1.0, x = 1.0, u;
int i = 0;
do {
i++;
x
*= 0.1;
u
= v + x;
cout << "i = "
<< i << " x = "
<<
setprecision(10) << x << " u = "
<<
setprecision(21) << u << endl;
}
while (u > v);
return
0;
}
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
上のプログラムをopenSUSE13.2で実行すると図１２の結果を得る。
図１２
２０桁ぐらいの精度であることがわかる。
long double型の値の範囲をdouble型のものと比較するプログラムを次のように作成した。
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
#include
<iostream>
#include
<iomanip>
#include
<cmath>
using
namespace std;
int
main(){
long
double v = 11000.0;
double
u = 600.0;
for (int i = 0; i
< 15; i++){
cout << "v = " << v <<
": " << setw(20) << left << exp(v)
<<
" u = " << u
<< ": " << exp(u) << endl;
v
+= 100.0;
u
+= 10.0;
}
return
0;
}
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
上のプログラムの実行結果を図１３に示す。
図１３
long double型のvは、double型のuより広い範囲の数値を表すことがわかる。
参考文献
岡本安晴（1998）「Delphiで学ぶデータ分析法」CQ出版
芝 祐順・南風原朝和（1990）「行動科学における統計解析法」東京大学出版会