【ブートストラップ法】標本平均の分散を評価してみる

背景

前回（【ブートストラップ法】経験分布関数について数値計算を通して理解を深める）は経験分布関数について数値計算で理解を深めました。今回はいよいよブートストラップ法を試してみることにします。

ブートストラップ法

ブートストラップ法は、未知の分布関数を経験分布関数で近似して統計的な推定を行う方法です。
経験分布関数からの乱数生成は、実は重複を許したリサンプリングになります（リサンプリングで作った標本をブートストラップ標本といいます）。これにより推定の安定性等を議論出来ます。

標本平均の分散を評価してみる

正規分布から$n$個の値を取り出し、母平均を推定することを考えます。この時問題になるのは標本平均の分散です。これについては$\sigma^2/n$となることを知っているので、母分散$\sigma^2$を不偏分散に置き換えれば推定することが出来ます。しかしもっと複雑な統計量については簡単にその分散を評価出来るとは限りません。やってみます。$\sigma^2=1$としました。

# include <random>
# include <iostream>

int main(int argc, char **argv)
{
	double m = 5;
	double s = 1;
	int n = 100;
	std::normal_distribution<> gauss(m,s);
	std::mt19937 mt(0);
	std::uniform_int_distribution<> u(0, n-1);

	std::vector<double> data(n);
	for (int i = 0; i < n; i++) data[i] = gauss(mt);
	double AV = 0;
	for (int i = 0; i < n; i++)
	{
		AV += data[i] / (double)n;
	}

	int B = atoi(argv[1]);//
	std::vector<double> avb(B);
	double av = 0;
	for (int b = 0; b < B; b++)
	{
		double sum = 0;
		for (int i = 0; i < n; i++)
		{
			double val = data[u(mt)];
			sum += val;
		}
		av += sum / (double)n;
		avb[b] = sum / (double)n;
	}
	av /= (double)B;
	double s2 = 0;
	for (int b = 0; b < B; b++)
	{
		s2 += (avb[b] - AV) * (avb[b] - AV);
	}
	std::cout << s2 / (double)(B-1) << std::endl;
	return 0;
}

結果を示します。

真の値は$\sigma^2/n=0.01$です。横軸は反復抽出の回数です。そこそこの値は出ていますが、ここから反復抽出の回数を増やしても真値に近づくわけではありません。経験分布関数を母分布と考えると、その分散$\hat\sigma^2$を標本サイズで割った値に近づきます。元のデータを増やせば経験分布関数が真の分布関数に近づくので精度がよくなります。

何が嬉しいのか

一般に統計量の期待値、分散の評価は手計算で簡単に出来るとは限りません。また、得られたデータの母分布の関数形が既知とも限りません。このような状況で母分布の関数形を仮定せず、手計算なしで統計量（推定量）の性質（不偏性や分散）の評価が出来そうです。