#はじめに

現代数理統計学の基礎 (共立講座数学の魅力)
2018/5/6、amazonで購入。
第６章「統計的推定」演習問題の問16に関するメモ
求まった最尤推定量が一致性を持つと記述されているが、実際は一致性を持たない。
自分の計算が本当に正しいかどうか不安だったので、Rプログラムを作成して検証してみた。

#問題の内容

$(X_1,Y_1),\dots,(X_n,Y_n)\ $を互いに独立な確率変数のペアとし、$X_i$と$Y_i$も独立でともに$N(\mu_i,\sigma^2)$に従うものとする。
(1) $\mu_i,\dots,\mu_n,\sigma^2$の最尤推定量を求めよ。
(2) $\sigma^2$の最尤推定量が一致性を持つことを示せ。
(3)【略】

#回答例

(1) 対数尤度を求めて最大化するパラメータを求める。

\begin{align}
f(x_i,y_i|\mu_i,\sigma^2)&=\frac{1}{2\pi\sigma^2}
\exp\left(-\frac{(x_i-\mu_i)^2+(y_i-\mu_i)^2}{2\sigma^2}\right)\\
\log\prod_{i=1}^n f(x_i,y_i)&=-n\log(2\pi)-n\log \sigma^2
-\frac{1}{2\sigma^2}\sum_{i=1}^n((x_i-\mu_i)^2+(y_i-\mu_i)^2)\\
\frac{\partial}{\partial \mu_i}\log\prod_{i=1}^n f(x_i,y_i)&=
\frac{1}{\sigma^2}(x_i-\mu_i+y_i-\mu_i)=0
\Rightarrow \mu_i=\frac{x_i+y_i}{2}\\
\therefore \quad T_i(X,Y)&=\frac{X_i+Y_i}{2}\quad{\mu_iの最尤推定量}\\
\log\prod_{i=1}^n f(x_i,y_i)&=-n\log(2\pi)-n\log \sigma^2
-\frac{1}{4\sigma^2}\sum_{i=1}^n(x_i-y_i)^2\\
\frac{\partial}{\partial \sigma^2}\log\prod_{i=1}^n f(x_i,y_i)&=
-\frac{n}{\sigma^2}+\frac{1}{4\sigma^4}\sum_{i=1}^n(x_i-y_i)^2=0
\Rightarrow \sigma^2=\frac{1}{4n}\sum_{i=1}^n(x_i-y_i)^2\\
\therefore \quad S(X,Y)&=\frac{1}{4n}\sum_{i=1}^n(X_i-Y_i)^2\quad{\sigma^2の最尤推定量}
\end{align}

演習問題の略解では、分母が4ではなく2となっている。何度も計算を確認したがやはり4だと思う。

(2) 得られた最尤推定量の収束値を求めてみる。

\begin{align}
X_i-Y_i&\sim N(0,2\sigma^2)\Rightarrow \frac{X_i-Y_i}{\sqrt{2}\sigma}\sim N(0,1)\\
\left(\frac{X_i-Y_i}{\sqrt{2}\sigma}\right)^2&\sim \chi^2_1
\Rightarrow E\left[\left(\frac{X_i-Y_i}{\sqrt{2}\sigma}\right)^2\right]=1\\
\therefore \quad S(X,Y)&\to_p\frac{\sigma^2}{2}(n\to\infty)\quad{大数の法則}
\end{align}

従って最尤推定量は一致性を持たない。真値の半分となっている。
期待値も$\frac{\sigma^2}{2}$であり、$n\to\infty$としても半分のままである。
分母が2なら一致性と不偏性を持つが、最尤推定量とは異なっている。

#プログラムによる検証

対数尤度を最大化するパラメータ$\sigma^2$が$\frac{1}{4n}\sum_{i=1}^n(x_i-y_i)^2$となることを確認する。

平均値$\mu_i$は一様乱数で生成。
optim関数に渡すパラメータは、$n$個の平均値$\mu_i$と$\sigma^2$から成る配列。
dnorm関数には標準偏差を指定する必要があるため、絶対値にしてから平方根とする。
optim関数は最小化するので、対数尤度の符号を変えることで最大化する。

set.seed(1)

sd = 2
n = 5
mu = runif(n)
xi = rnorm(n,mean=mu,sd=sd)
yi = rnorm(n,mean=mu,sd=sd)
mu_mle = (xi+yi)/2
s2_mle = sum((xi-yi)^2)/4/n

my_logLik = function(pa){
  mn = pa[1:n]
  s1 = sqrt(abs(pa[n+1]))
  -sum(dnorm(x=xi,mean=mn,sd=s1,log=T)+dnorm(x=yi,mean=mn,sd=s1,log=T))
}

par = rep(1,n+1)
rt = optim(par,my_logLik,control=list(maxit=10000))

結果は以下の通り。

> c(s2_mle,rt$par[n+1])
[1] 2.328907 2.328533
> 
> data.frame(mu_mle,par=rt$par[1:n])
      mu_mle        par
1  1.5321708  1.5321919
2  3.1914187  3.1910387
3 -0.2035032 -0.2028649
4 -0.8193685 -0.8192448
5 -1.2406955 -1.2406410

対数尤度を最大化するパラメータとして、$\sigma^2$は$\frac{1}{4n}\sum_{i=1}^n(x_i-y_i)^2$、$\mu_i$は$\frac{x_i+y_i}{2}$とほぼ一致する値が得られた。

#おわりに

【注意事項】本メモの内容が本当に正しいかどうかは分かりません。