【区間推定】信頼区間の解釈、構成法について数値実験してみる

背景

統計学についてつまみ食いの知識しかないのはいけないと思い、「現代数理統計学」(竹村彰通)（以後教科書）を読み進めています。１週目は興味のある所を重点的に読んでいます。昨日から第9章の区間推定に入りました。
信頼区間の解釈については世間でも誤解が多く、私自身もよく理解出来ていませんでした（そもそも普段区間推定を使わない）。そこで、信頼区間の解釈、構成法について数値実験を通して理解を深めることにしました。

信頼区間は統計量

信頼区間の定義は教科書によれば

P_\theta[L(X) < \theta < U(X)] \geq 1-\alpha

という性質を持つ区間として定義されています。$L,R$はそれぞれLower、Upperに由来すると思われます。
$L(X),R(X)$と書かれている通り、これらは統計量、つまり標本の関数、確率的に変動する値です。データのサンプリングをし、統計量$L(X),R(X)$の値を算出することを無限に繰り返します。このとき$1-\alpha$の確率で真の母数$\theta$を$L,R$の間に含むことを上の式は意味します。
つまり信頼区間とは言っても範囲内に必ずしも真の母数が含まれているとは限りません。例えば$\alpha=0.05$であれば、試行を100回繰り返した場合（つまり計算される信頼区間も100個ある）、95回くらいはその信頼区間に真の母数を含みますが、5回くらいは真の母数が外れてしまっていることになります。

信頼区間の構成法

では、どのようにして信頼区間を構成すればよいのでしょう。実は両側検定を利用することで信頼区間を構成出来ます。帰無仮説は$\theta=\theta_0$、対立仮説は$\theta\neq\theta_0$です。今、標本$X$があるとします。これはある帰無仮説$\theta=\theta_0$では受容されるかもしれませんが、別の帰無仮説$\theta=\theta_0'$では棄却されるかもしれません。与えられた$X$が有意水準$\alpha$で受容される$\theta$の集合を考えると、実はそれが信頼区間となっています。
どのような検定を選ぶかによって算出される信頼区間は異なります。

数値実験

ポアソン分布から$n$個の値をサンプリングして$\lambda$の区間推定をしてみます。今回は簡単に尤度比検定を利用して信頼区間を構成してみます。前回の記事（【検定論】尤度比検定を数値実験してみる）によれば対数尤度比は

\log L = n(\lambda_0-\bar X) + \sum_i X_i\log\frac{\bar X}{\lambda_0}

であり、有意水準5%の検定は

2\log L \leq \chi_{0.05}^2(1)\simeq 3.84 \to {\rm accept}

で与えられます。上の式で$\lambda_0$の値を振って不等式を満たす範囲を算出すればそれが信頼区間です。このような算出方法で本当に信頼区間の条件を満たすのでしょうか。$\lambda_0$の値を真の値に設定し（つまり帰無仮説が正しいとし）、$X$のサンプリングを繰り返すことを考えます。このとき、有意水準5%の検定であれば確率95%で受容されます。つまり、信頼区間の算出を繰り返すと95%は真の母数を含むことになり、信頼区間の定義を満たします。
今回は標本サイズ$n=100,\lambda=3$として数値実験します。$\lambda_0$の値を$0.01$ずつ増やして信頼区間を算出してみます。これを$100$回繰り返します。

#include <iostream>
#include <fstream>
#include <random>
#include <vector>

int factorial(int k)
{
	int ret = 1;
	for (int i = 1; i <= k; i++)
	{
		ret *= i;
	}
	return ret;
}

int main(int argc, char **argv)
{
	std::mt19937 mt(0);
	double lambda = 3;
	int N = 100;//標本サイズ
	std::poisson_distribution<int> poi(lambda);
	double alpha = 0.05;
	double chi2 = 3.84;
	
	int Iter = 100;
	int reject = 0;
	for (int iter = 0; iter < Iter; iter++)
	{
		std::vector<int> data(N);
		double av = 0;
		for (int n = 0; n < N; n++)
		{
			data[n] = poi(mt);
			av += data[n] / (double)N;
		}
		double L, R;
		for (double lambda0 = 0; lambda0 < 20; lambda0 += 0.01)
		{
			double ll = N * (lambda0 - av + av * log(av / lambda0));
			if (2 * ll <= chi2)
			{
				L = lambda0;
				break;
			}

		}
		for (double lambda0 = 20; lambda0 >=0; lambda0 -= 0.01)
		{
			double ll = N * (lambda0 - av + av * log(av / lambda0));
			if (2 * ll <= chi2)
			{
				R = lambda0;
				break;
			}
		}
		std::cout << iter << " " << av << " " << L << " " << R << std::endl;
	}

	return 0;
}

大半は$\lambda=3$を含んでいますが、いくつか外れてしまっているものもありますね。ちなみに今回は尤度比検定を行ったので、正確に95%信頼区間とはなっていません。
ちなみに標本サイズを$n=1000$まで増やすと信頼区間の幅がぐっと小さくなります。

感想

信頼区間とは何かがわかってきました。目の前の限られたデータセットから95%信頼区間を構成したとして、運悪く残りの5%に当たってしまっていたらショックですね。その場合は97%信頼区間等にすればよいのでしょうが、当然幅は広がっていくわけで。もどかしい。