背景
前回（【検定論】片側検定問題の一様最強力検定）は単純検定のみならず、片側検定においてもある条件を満たせば一様最強力検定を構成出来ることを述べました。
しかし一般には一様最強力検定が存在するという保証はありません。今回数値実験する尤度比検定は理論的保証があるわけではありませんが、合理的な検定が得られる場合が多く、実用的なのでよく用いられるそうです。（参考：「現代数理統計学」竹村彰通）
どのように検定すればよいか
帰無仮説を$\theta\in\Theta_0$、対立仮説を$\theta\in\Theta_1$とします。尤度比は
L = \frac{\max_{\theta\in \Theta_1}f(x,\theta)}{\max_{\theta\in \Theta_0}f(x,\theta)}
で定義されます。ネイマン・ピアソンの補題同様、ある$c$を考え、$L>c$のとき棄却します。$c$は有意水準$\alpha$から計算されますが、一般に簡単に求まるとは限りません。しかし、帰無仮説の下で自由に動かせるパラメータの数を$q$、対立仮説の下で自由に動かせるパラメータの数を$p+q$とすると、標本サイズ$n\to\infty$のとき
2\log L > \chi_\alpha^2(p) \to {\rm reject}
が有意水準$\alpha$の検定となります。これなら機械的にできますね。
数値実験
ポアソン分布から$n$個の値をサンプリングすることを考えます。帰無仮説は$\lambda=\lambda_0$、対立仮説は$\lambda\neq\lambda_0$とします。有意水準は5%とします。まず、帰無仮説の下で尤度を最大とするのは言わずもがな$\hat\lambda=\lambda_0$です。対立仮説の下で尤度を最大とするのは$\hat\lambda=\bar {X}$です。よって
L = \frac{\prod_i e^{-\bar X}\bar X^{x_i}}{\prod_i e^{-\lambda_0}\lambda_0^{x_i}}
となります。対数尤度比は
\log L = n(\lambda_0-\bar X) + \sum_i X_i\log\frac{\bar X}{\lambda_0}
となります。帰無仮説で動かせるパラメータは$q=0$、対立仮説で自由に動かせるパラメータは$p+q=1$ですので、尤度比検定は
2\log L > \chi_{0.05}^2(1)\simeq 3.84 \to {\rm reject}
となります。
それでは$\lambda_0=3$として数値実験してみます。
#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 lambda0 = 3;
double lambda = 3;
int N = atoi(argv[1]);//標本サイズ
std::poisson_distribution<int> poi(lambda);
double alpha = 0.05;
double chi2 = 3.84;
int Iter = 10000;
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 ll = N * (lambda0-av);
for (int n = 0; n < N; n++)
{
ll += data[n] * log(av / lambda0);
}
if (2*ll > chi2) reject++;
}
std::cout << reject / (double)Iter << std::endl;
return 0;
}
標本サイズ10以下は怪しいですが、それ以上では有意水準5%前後の値となっています。
感想
とりあえず困ったら尤度比とっとけ、って感じでしょうか。