배경

「현대수리 통계학」(다케무라 아키히로)(이후, 교과서)는 은근하게 흥미있는 곳은 읽었습니다. 제11장 제2절의, 회귀 모델을 충분히 통계량의 관점에서 고찰하고 있는 곳이 재미있었습니다.
2주째는 세세한 증명도 쫓을 예정입니다만, 그와 동시에 다른 흥미있는 화제도 쫓고 싶다고 하는 것으로 부트스트랩법입니다. 구간 추정도 할 수 있으므로, 교과서에서 공부한 방법과의 비교도 하고 싶습니다만, 이번은 준비로서 경험 분포 함수에 대해서, 수계산과 수치 계산으로 이해를 깊게 합니다.

경험 분포 함수란?

누적 도수 분포라고 하는 것이 익숙해질지도 모릅니다. 경험 분포 함수는 표본에 의존하기 때문에 통계량이며 누적 분포 함수의 추정량입니다. 정의는

\hat F_x(\{x\}) = \frac{1}{n}\sum_i I(x_i\leq x)

입니다. 단 $I$는 지시 함수로

I(x_i\leq x) = 
  \left\{
    \begin{array}{ll}
      1 & x_i \leq x \\
      0 & {\rm otherwise}
    \end{array}
  \right.

입니다.
포아송 분포로 해보자. $\lambda=5$로 설정했습니다.

표본 크기를 늘리면 진정한 분포 함수에 가깝다는 것을 알 수 있습니다. 코드는 나중에 표시됩니다.

기대치, 분산을 계산해 본다

그렇다면 얼마나 많은 샘플 크기가 필요합니까? 경험 분포 함수는 통계량이므로 기대치, 분산을 계산할 수 있습니다. 해보자. 이산 분포로 생각합니다.

\begin{eqnarray}
{\rm E}[\hat F_\lambda(\{x\})] &=& \sum_{\{x\}}\hat F_\lambda(\{x\})p(\{x\})\\
&=& \sum_{x_0}\cdots\sum_{x_{n-1}} \hat F_\lambda(\{x\})p(\{x\})\\
&=& \sum_{x_0}\cdots\sum_{x_{n-1}} \frac{1}{n}\sum_{i=0}^{n-1}I(x_i\leq\lambda)p(x_0)\cdots p(x_{n-1})\\
&=& \frac{1}{n} \sum_{x_0}\cdots\sum_{x_{n-1}}\{I(x_0\leq\lambda)+\cdots\}p(x_0)\cdots p(x_{n-1})
\end{eqnarray}

되고,

\sum_{x_0}\cdots\sum_{x_{n-1}}I(x_0\leq\lambda)p(x_0)\cdots p(x_{n-1}) = \sum_{x_0\leq\lambda}p(x_0) = F(\lambda)

보다

{\rm E}[\hat F_\lambda(\{x\})] = F(\lambda)

됩니다. 다음으로 분산을 생각합니다.

\begin{eqnarray}{\rm V}[\hat F_\lambda(\{x\})] &=& {\rm E}[\hat F_\lambda^2] - {\rm E}[\hat F_\lambda]^2\end{eqnarray}

보다 제1항을 생각합니다.

\begin{eqnarray}{\rm E}[\hat F_\lambda^2] &=& \frac{1}{n^2}\sum_{x_0x_1\cdots}\sum_{ij}I(x_i\leq\lambda)I(x_j\leq\lambda)p(x_0)\cdots p(x_{n-1})\end{eqnarray}

라고 쓸 수 있습니다.

\begin{eqnarray}\sum_{x_0x_1\cdots}I(x_i\leq\lambda)I(x_j\leq\lambda)p(x_0)\cdots p(x_{n-1}) &=& \left\{    \begin{array}{ll}      F(\lambda) & i=j \\      F(\lambda)^2 & i\neq j    \end{array}  \right.\end{eqnarray}

이고, $i=j$가 되는 경우가 $n$대로, $n\neq j$가 되는 경우가 $n(n-1)$인 것에 주의하면,

{\rm V}[\hat F_\lambda(\{x\})] = F(\lambda)(1-F(\lambda))/n

됩니다. $n\to\infty$에서 0이 될 수 있습니다. 또, $F=0.5$일 때가 가장 큰 값을 취하는 것을 알 수 있습니다.
분산의 $n$ 의존성을 수치 계산으로 확인해 보겠습니다. 어머니 평균 $5$, $F(5)$의 추정을 실시했습니다.

깔끔하게 일치합니다.
이번에 사용한 코드는 다음과 같습니다.

#include <random>
#include <iostream>
#include <vector>
#include <algorithm>

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 m = 5;//母数

    std::poisson_distribution<> poi(m);

    int Iter = 100000;
    std::vector<double> est;
    int lambda0 = m;
    int n = atoi(argv[1]);

    auto f = [&](int lambda)
    {
        return exp(-m) * pow(m, lambda) / factorial(lambda);
    };
    auto F = [&](int lambda)//分布関数
    {
        double ret = 0;
        for (int X = 0; X <= lambda; X++)
        {
            ret += f(X);
        }
        return ret;
    };

    for (int iter = 0; iter < Iter; iter++)
    {

        std::vector<int> x(n);
        for (int i = 0; i < n; i++)
        {
            x[i] = poi(mt);
        }
        std::sort(x.begin(), x.end());

        auto Fe = [&](int lambda)//経験分布関数
        {
            double ret = 0;
            for (const auto& e : x)
            {
                if (e > lambda) break;
                ret += 1. / (double)n;
            }
            return ret;
        };

        est.push_back(Fe(lambda0));
        double sum = 0;
        /*for (int X = 0;  X <= 2 * m; X++)
         *      {
         *                  sum += exp(-m) * pow(m, X) / factorial(X);
         *                              std::cout << X << " " << Fe(X) << " " << F(X) << std::endl;
         *                                      }*/
    }

    double estav = 0,s2 = 0;
    for (int iter = 0; iter < Iter; iter++)
    {
        estav += est[iter] / (double)Iter;
    }
    for (int iter = 0; iter < Iter; iter++)
    {
        s2 += (est[iter] - estav)* (est[iter] - estav) / (double)Iter;
    }

    std::cout << n << " " << s2 << " " << F(lambda0) * (1. - F(lambda0)) / (double)n << std::endl;;
    //std::cout << estav << " " << f(lambda0);
    //  
    return 0;
}

향후 예정

신뢰 구간의 추정 등을 하고 싶습니다.