稲垣『数理統計学』演習問題8.3

最近数理統計を復習している。というか理解があまりにあやふやなので普通に勉強している。
そういう流れで記事タイトルの問題をゴリ押しの解法で解いたのだが、あまりに面倒だった。計算過程を忘れる前にメモっておく。要領のいい別解があるかもしれないのでそのようなコメントも期待しつつ。

問題8.3
成功の確率が \displaystyle p \displaystyle n個の独立なベルヌーイ試行 \displaystyle  \varepsilon _1 , \varepsilon _2 , \dots , \varepsilon _nによる母分散  \displaystyle p(1-p) の推定量 \displaystyle T = c\overline{X} (1 - \overline{X}) とする。ここで、
 \displaystyle \overline{X} = \frac{1}{n}(\varepsilon _1 + \varepsilon _2 + \cdots + \varepsilon _n):標本比率である。  \displaystyle Tが母分散の不偏推定量であるとき、 \displaystyle cの値を求めよ。さらに、そのときの推定量の分散を求めよ。

解答
 \displaystyle Tが母分散の不偏推定量であるには、 \displaystyle  E[T] =p(1-p) であればよい。
 \displaystyle  E[T] = cE[\overline{X}] - cE[\overline{X} ^2] である。 \displaystyle  E[\overline{X}] = \sum_{i=1}^n \frac{1}{n} E[\varepsilon _i ] = p で、
 \displaystyle  E[\overline{X} ^2 ] = \sum_{i,j=1}^n \frac{1}{n ^2} E[\varepsilon _i \varepsilon _j ] であるが、これら \displaystyle  n^2 項のうち、
   \displaystyle  \varepsilon _i ^2のパターン…  \displaystyle n個で、期待値 \displaystyle  p
   \displaystyle  \varepsilon _i \varepsilon _j (i \not = j)のパターン…  \displaystyle n(n-1)個で、期待値 \displaystyle  p ^2
なので \displaystyle  E[\overline{X} ^2 ] =\frac{1}{n ^2} (n p + (n ^2 - n) p ^2)
代入して \displaystyle  c について解けば、不偏推定量になるのは  \displaystyle  c = \frac{n}{n-1}

このときの分散を求める。そのために、まず \displaystyle  E[T ^2 ] を求めよう。
 \displaystyle  E[T ^2 ] = c ^2 E[\overline{X} ^2 (1 - \overline{X} ) ^2 ] = c ^2(E[ \overline{X} ^2 ] - 2 E[ \overline{X} ^3 ] +  E[ \overline{X} ^4 ]) である。
先程求めたように、 \displaystyle  E[\overline{X} ^2 ] =\frac{1}{n ^2} (n p + (n ^2 - n) p ^2)
 \displaystyle  E[\overline{X} ^3 ] = \sum_{i,j,k=1}^n \frac{1}{n ^3} E[\varepsilon _i \varepsilon _j \varepsilon _k ] であるが、これら \displaystyle  n^3 項のうち、
   \displaystyle  \varepsilon _i ^3のパターン…  \displaystyle n個で、期待値 \displaystyle  p
   \displaystyle  \varepsilon _i ^2 \varepsilon _jのパターン…  \displaystyle 3n(n-1)個で、期待値 \displaystyle  p ^2
   \displaystyle  \varepsilon _i \varepsilon _j \varepsilon _k のパターン…  \displaystyle n(n-1)(n-2)個で、期待値 \displaystyle  p ^3
より \displaystyle  E[\overline{X} ^3 ] =\frac{1}{n ^3} (n p + 3 n (n - 1) p ^2 + n (n - 1) (n - 2) p ^3 )
 \displaystyle  E[\overline{X} ^4 ] = \sum_{i,j,k,l=1}^n \frac{1}{n ^4} E[\varepsilon _i \varepsilon _j \varepsilon _k \varepsilon _l ] であるが、これら \displaystyle  n^4 項のうち、
   \displaystyle  \varepsilon _i ^4のパターン…  \displaystyle n個で、期待値 \displaystyle  p
   \displaystyle  \varepsilon _i ^3 \varepsilon _jのパターン…  \displaystyle 4n(n-1)個で、期待値 \displaystyle  p ^2
   \displaystyle  \varepsilon _i ^2 \varepsilon _j ^2のパターン…  \displaystyle 3 n(n-1)個で(ダブルカウントに注意!!)、期待値 \displaystyle  p ^2
   \displaystyle  \varepsilon _i ^2 \varepsilon _j \varepsilon _k のパターン…  \displaystyle 6n(n-1)(n-2)個で(同上)、期待値 \displaystyle  p ^3
   \displaystyle  \varepsilon _i \varepsilon _j \varepsilon _k \varepsilon _l のパターン…  \displaystyle n(n-1)(n-2)(n-3)個で、期待値 \displaystyle  p ^4
より
 \displaystyle  E[\overline{X} ^4 ] =\frac{1}{n ^4} (n p + 7 n (n - 1) p ^2 + 6 n (n - 1 ) (n - 2) p ^3 + n (n - 1) (n -2 ) (n - 3) p ^4)
と求まる。これらを代入して整理すると(この一言にあまりにもつらい計算が込められている)、
 \displaystyle  E[T ^2 ] = \frac{p (p - 1)}{n (n - 1)} \{ (n ^2 - 5 n + 6) (p ^2 - p) - n + 1 \}
よって、
  \displaystyle  \begin{array}\displaystyle  V[ T] &= E[T ^2 ] - E[T ] ^2 = E[T ^2 ] - p ^2 (1 - p) ^2 \\
&=\displaystyle    \frac{p (p - 1)}{n (n - 1)} \{ ( - 4 n + 6) (p ^2 - p) - n + 1 \} \\
&=\displaystyle    \frac{p (p - 1)}{n (n - 1)} \{ - n (4 p ^2 - 4 p + 1) + (4 p ^2 - 4 p + 1) + 2 p ^2 - 2 p \} \\
&=\displaystyle    \frac{p (p - 1)}{n (n - 1)} \{ - (n - 1) (2 p - 1) ^2 + 2 p (p - 1) \} \\
&=\displaystyle    \frac{p (1 - p)}{n} \left\{ (2 p - 1) ^2 + 2 \frac{p (1 - p)}{n - 1}  \right\}
\end{array}
と求まる。■

一様乱数の和が1を超える個数の期待値

はてなブログtex記法を思い出すためにひとつ書こうと思う。140字以上の文章をインターネットに載せるのは久しぶりだ。


twitterで流れてきた問題の解答。とりうる確率変数の値が無限個であるが、確率変数がそれぞれの値をとる確率を求めて、無限級数の形で期待値を計算するだけの愚直なアイデアで解ける。途中で登場する積分計算が難所だろう。要はn次元単体の体積なのだが、実際に計算するためには一捻り要る。

[解答]
はじめて和が1をこえる項数が nである確率をP_nとおく。求めたい期待値は \sum_{n=1}^\infty nP_nと表される。
\triangle_n(t) := \{ (x_1,x_2,\dots , x_n) \in \mathbb{R}^n \mid 0 \leq x_i \leq t, i=1,2,\dots ,n, 0 \leq x_1+x_2+\dots +x_n \leq t \}
とおく。また、 Q_n := P(a_1+a_2+ \cdots +a_n \leq 1)とおく。
 a_iは一様分布U(0,1)にしたがうので、正の整数 nに対して、
 \begin{array}. Q_{n+1} &= \int_{\triangle_{n+1}(1)}dx_1dx_2\cdots dx_{n+1} \\ 
&= \int_0^1 \left( \int_{\triangle_n(1-x_{n+1})}dx_1dx_2 \cdots dx_n \right) dx_{n+1}\ (\text{Fubiniの定理より}) \\
&= \int_0^1 \left( \int_{\triangle_n(1)}(1-x_{n+1})^n dx_1dx_2 \cdots dx_n \right) dx_{n+1} \\
&= \int_0^1 (1-x_{n+1})^n Q_n dx_{n+1} \\
&= \frac{1}{n+1}Q_n \end{array}
これと Q_1 = 1から帰納的に Q_n = \frac{1}{n!}だとわかる。Q_nははじめて和が1をこえる項数がnより大きい確率なので、 n \geq 2に対して
 P_n = Q_{n-1} - Q{n} = \frac{1}{(n-1)!} - \frac{1}{n!}
これはn=1のときにも成り立つので、求める期待値は
 \begin{array}
.\sum_{n=1}^\infty nP_n &= \sum_{n=1}^\infty n \left( \frac{1}{(n-1)!} - \frac{1}{n!} \right) \\
&= \sum_{n=1}^\infty \left( \frac{n-1}{(n-1)!} + \frac{1}{(n-1)!} - \frac{n}{n!} \right) \\
&= \sum_{n=2}^\infty \frac{1}{(n-2)!} \\
&= e.\ ■\end{array}

論より証拠、実際に計算してみよう(計算機上でだが)。

set.seed(114514)

NUM <- 10000000

trial <- function(){
  count <- 0
  S <- 0
  while(S <= 1){
    count <- count + 1
    X <- runif(1,0,1)
    S <- S + X
  }
  return(count)
}

E <- 0
for(i in 1:NUM){
  E <- E + trial()
}
E <- E / NUM

print(E)

実行結果は以下:

[1] 2.718196

ちょっとだけeに近いが、あまり精度はよくない。