背景
「現代数理統計学」(竹村彰通)（以後、教科書）はさらっと興味ある所は読みました。第１１章第２節の、回帰モデルを十分統計量の観点から考察しているところが面白かったです。
２週目は細かい証明も追う予定ですが、それと同時に他の興味ある話題も追いたい、ということでブートストラップ法です。区間推定も出来るので、教科書で勉強した方法との比較もしたいと思っていますが、今回は準備として経験分布関数について、手計算と数値計算で理解を深めます。
経験分布関数とは
累積度数分布と言った方が馴染みがあるかもしれません。経験分布関数は標本に依存しますので、統計量であり、累積分布関数の推定量です。定義は
\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}
より第一項を考えます。
\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$でゼロになることがわかります。また、$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;
}
今後の予定
信頼区間の推定などをしたいです。