背景

検定論でネイマン・ピアソンの補題というものが出てきます。これは尤度比検定が優れた検定であることを保証するものです。式を追ってるだけだとピンとこないので数値実験してみます。

ネイマン・ピアソンの補題

帰無仮説は母数が$\theta=\theta_0$、対立仮説は$\theta=\theta_1$とします。尤度比
$$
f(x,\theta_1)/f(x,\theta_0)
$$
の値を計算し、ある値$c$よりも大きいとき帰無仮説を棄却します。$c$は有意水準$\alpha$の値によって設定できます。連続分布の場合、帰無仮説が正しいときこれを棄却する確率は$c$を調整することで$\alpha$に一致させられるので、これがそのまま最強力検定となります。離散分布の場合は一般に一致させられないため、尤度比が$c$に一致するときはコインを投げて受容するか棄却するかを決めます。

数値実験

ポアソン分布から値をひとつだけサンプリングすることを考えます。帰無仮説は$\lambda=\lambda_0$、対立仮説は$\lambda=\lambda_1\;(\lambda_1>\lambda_0)$とします。尤度比は
$$
\frac{e^{-\lambda_1}\lambda_1^x}{e^{-\lambda_0}\lambda_0^x}
$$
となります。
$$
\frac{e^{-\lambda_1}\lambda_1^x}{e^{-\lambda_0}\lambda_0^x} > c
$$
のとき帰無仮説を棄却します。両辺の対数を取ってやれば結局
$$
x > c'
$$
のとき棄却すればよいことがわかります。有意水準を$\alpha$以下にしたいため、帰無仮説が正しいとき
$$
P_{\lambda_0}(x>c') \leq \alpha
$$
でなければなりません。等号が成立すれば検定にコイントスを用いる必要はないのですが、
$$
\alpha - P_{\lambda_0}(x>c') \equiv \beta
$$
が$0$でない場合はコイントスを導入する必要があります。これは帰無仮説を必要以上に受容しているためです。これを調整するため、尤度比が$c'$に一致する場合は確率$r$で帰無仮説を棄却します。$r$は
$$
P_{\lambda_0}(x=c') \times r = \beta
$$
から決定できます。
それでは具体的な値を用いて検定してみます。$\lambda_0=3,\lambda_1=5$とします。有意水準は5%以下とします。最初に$c'$の値を決めます。

#include <iostream>
#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 N = 10;
    std::vector<double> p(N);//帰無仮説の下での母分布
    double lambda0 = 3;
    double sum = 0;
    int c = 6;
    double alpha = 0.05, beta;
    for (int k = 0; k <N; k++)
    {
        p[k] = exp(-lambda0) * pow(lambda0, k) / factorial(k);
        sum += p[k];
        std::cout << "x > " << k << " : " << 1-sum << std::endl;
    }
    return 0;
}

x > 0 : 0.950213
x > 1 : 0.800852
x > 2 : 0.57681
x > 3 : 0.352768
x > 4 : 0.184737
x > 5 : 0.0839179
x > 6 : 0.0335085
x > 7 : 0.0119045
x > 8 : 0.00380299
x > 9 : 0.00110249

これから$c'$は6以上でなければならないことがわかります。$c'=6$とすれば$r\simeq0.33$となります。これ以外の値に設定すると$r$は1を超え、確率化検定(コイントスを導入した検定)が出来なくなり、有意水準は満たしますが、最強力検定となりません。一般に$c$の値は一意に決定できます。では検定をしてみます。

#include <iostream>
#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 N = 10;
    std::vector<double> p(N);//帰無仮説の下での母分布
    double lambda0 = 3;
    double sum = 0;
    int c = 6;
    double alpha = 0.05, beta;
    for (int k = 0; k <N; k++)
    {
        p[k] = exp(-lambda0) * pow(lambda0, k) / factorial(k);
        sum += p[k];
        std::cout << "x > " << k << " : " << 1-sum << std::endl;
        if (k  == c) beta = alpha - (1-sum);
    }
    double r = beta / p[c];
    std::cout <<"r = "  << r << std::endl;

    int Iter = 1000000;
    std::uniform_real_distribution<> u(0, 1);
    std::mt19937 mt(0);

    int reject = 0;
    double lambda1 = 5;
    std::poisson_distribution<> poi(lambda0);//検出力を評価するときはlambda0 -> lambda1
    for (int iter = 0; iter < Iter; iter++)
    {
        int x = poi(mt);

        if (x > c) reject++;
        else if (x == c)
        {
            if (u(mt) < r) reject++;//確率化検定をしないときはコメントアウト
        }
    }
    std::cout << reject / (double)Iter;
    return 0;
}

有意水準：0.050166

有意水準5%を概ね満たしていることが確認できました。検出力は約$0.286$です。確率化検定をしない場合は

有意水準：0.03364

であり、有意水準5%以下は満たしていますが、検出力は約$0.238$で確率化検定よりも低い値となっていることが確認できます。

まとめ

ネイマン・ピアソンの補題を数値実験で確認しました。確率化決定関数なんていつ使うんだ、って思ってましたがこんなとこで出会うとは。
学生時代、検定論は試験をパスするために少し勉強しただけでしたが、勉強してみると結構面白いですね。