最尤法
確率密度関数 $f(x; \theta)$ から独立に $n$ 個の標本 $x_1,...,x_n$ を得たとする. 標本は独立であることから, 同時確率関数を $L(\theta)$ とすると以下のようになる.
\begin{align}
L(\theta) &= f(x_1,...x_n; \theta) \\
&= f(x_1; \theta) \times ... \times f(x_n; \theta) \\
&= \prod_{i=1}^{n} f(x_i; \theta)
\end{align}
$L(\theta)$ のうち標本 $x_1,...,x_n$ は観測済みの固定された値であるので, $L(\theta)$ は $\theta$ の関数である.
最尤法とは, この尤度関数 $L(\theta)$ が最大となるような $\theta$ の値を求める手法である. つまり, 得られた標本が最も「起こりやすい」ようなパラメータ $\theta$ を選ぶ方法である. このときの $\theta$ の値を 最尤推定量 と呼ぶ.
なお, 尤度関数 $L(\theta)$ の対数を取ることで計算がしやすくなることが多く,
尤度関数 $L(\theta)$ の対数を取った $l(\theta)$ は対数尤度と呼ばれる.
\begin{align}
l(\theta) &= \log L(\theta) \\
&= \log \prod_{i=1}^{n} f(x_i; \theta) \\
&= \sum_{i=1}^{n} \log f(x_i; \theta)
\end{align}
例2
対数尤度を計算すると,
\begin{align}
l(\mu, v) &= \log \prod_{i=1}^{n} \frac {1} {\sqrt{2\pi v}} \exp \bigg\{ - \frac {(x_i - \mu)^2} {2v} \bigg\} \\
&= \sum_{i=1}^{n} \log \frac {1} {\sqrt{2\pi v}} \exp \bigg\{ - \frac {(x_i - \mu)^2} {2v} \bigg\} \\
&= \sum_{i=1}^{n} \Bigg\{ \log \frac {1} {\sqrt{2\pi v}} + \bigg\{ - \frac {(x_i - \mu)^2} {2v} \bigg\} \Bigg\} \\
&= - \frac {n} {2} \log (2\pi v) - \frac {1} {2v} \sum_{i=1}^{n} (x_i - \mu)^2 \\
\end{align}
ワークブックとは求め方が違うが, $\partial l(\mu, v) / \partial \mu=0$ となる $\mu$ が $l(\mu, v)$ を最大にするので,
\begin{align}
\frac {\partial l(\mu, v)} {\partial \mu} &= \frac {1} {v} \sum_{i=1}^{n} (x_i - \mu) = 0 \\
\mu &= \frac {1} {n} \sum_{i=1}^{n} x_i = \bar{x}
\end{align}
同様に, $v$ についても
\begin{align}
\frac {\partial l(\mu, v)} {\partial v} &= - \frac {n} {2 v} + \frac {\sum_{i=1}^{n} (x_i - \mu)^2} {2v^2} = 0 \\
\mu &= \frac {1} {n} \sum_{i=1}^{n} (x_i - \mu)^2 = s
\end{align}
となり, $\mu, v$ は標本平均 $\bar{x}$, 標本分散 $s$ のとき最尤推定値となる.
モーメント法
$m$ 個のパラメータ $\theta = (\theta_1, ..., \theta_m)$ をもつ確率密度関数 $f(x; \theta)$ を考える.
$1,2,...,k$ 次の中心モーメントは定義より以下のようになる.
$$
\mu_1 = \int x f(x; \theta) dx
$$
$$
\mu_2 = \int (x - \mu_1)^2 f(x; \theta) dx
$$
$$
\mu_k = \int (x - \mu_1)^k f(x; \theta) dx
$$
これらの各左辺はパラメータ $\theta$ の関数であるのでこれを $m_k(\theta)$ と書くと, $\mu_k = m_k(\theta)$ のように表される. これらは $\theta$ についての連立方程式であるので, 適当な関数 $g_j$ を用いて
$$
\theta_j = g_j(\mu_1, ..., \mu_k), j=1,...,m
$$
のように $\theta$ を求めることができる.
モーメント法では, $\mu_1$ が期待値, $\mu_2$ が分散であるから, 標本平均や標本分散のような推定量を, $\mu_k$ の推定量とすることで, $\theta$ を求める.
つまり,
$$
\mu_1 = \hat{\mu_1} = \bar{X} = \frac {1} {n} \sum_{i=1}^{n} X_i
$$
$$
\mu_2 = \hat{\mu_2} = \frac {1} {n} \sum_{i=1}^{n} (X_i - \bar{X})^2
$$
$$
\mu_k = \hat{\mu_k} = \frac {1} {n} \sum_{i=1}^{n} (X_i - \bar{X})^k
$$
として,
$$
\hat{\theta_j} = g_j(\hat{\mu_1}, ..., \hat{\mu_k}), j=1,...,m
$$
のように $\theta$ を求める
例3
$Ga(\alpha, 1/\lambda)$ の平均と分散は
$$
\mu_1 = \alpha / \lambda
$$,
$$
\mu_2 = \alpha / \lambda^2
$$
である.
$\mu_1, \mu_2$ の推定量として標本平均 $\bar{X}$, 標本分散 $S$ を用いると,
$$
\bar{X} = \alpha / \lambda
$$,
$$
S = \alpha / \lambda^2
$$
これを$\alpha, \lambda$ について解くと, モーメント法による推定量 $\alpha = \mu_1^2/\mu_2, \lambda=\mu_1/\mu_2$ が得られる.
推定量の評価指標
平均二乗誤差
$$
MSE_{\theta}(\hat{\theta}) = E_{\theta}[(\hat{\theta} - \theta)^2]
$$
相対効率 : 二つの推定量の推定精度の比較, $\hat{\theta}_1$ の $\hat{\theta}$ に対する相対効率
$$
e(\hat{\theta}_1, \hat{\theta}_2) = \frac {E[(\hat{\theta}_1 - \theta)^2]} {E[(\hat{\theta}_2 - \theta)^2]}
$$
点推定の性質
不偏推定量
$E[\hat{\theta}] = \theta$ となる推定量を不偏推定量という.
$b(\theta) = E[\hat{\theta}] - \theta$ を推定量 $\hat{\theta}$ のバイアスと呼ぶので, 不偏推定量はバイアスが0の推定量と言い換えられる.
平均二乗誤差 $E[(\hat{\theta} - \theta)^2]$ は $E[X^2] = E[X]^2 + V[X]$ と真のパラメータ $\theta$ は定数であることから以下のように分けることができる.
\begin{align}
E[(\hat{\theta} - \theta)^2] &= (E[\hat{\theta} - \theta])^2 + V[\hat{\theta} - \theta] \\
&= (E[\hat{\theta}] - \theta)^2 + V[\hat{\theta}] \\
&= b(\theta) + V[\hat{\theta}]
\end{align}
これを平均二乗誤差のバイアス・バリアンス分解という. バイアス項が0となる不偏推定量に限ったうえで, 平均二乗誤差を最小化する推定量は, 分散 $ V[\hat{\theta}]$ を最小化する推定量にほかならない. このような推定量を, 一様最小分散不偏推定量という.
有効推定量
以下, クラメールラオの不等式の等号を満たすような不偏推定量を有効推定量という. よって, 有効推定量ならば一様最小分散不偏推定量である.
クラメールラオの不等式の前に, フィッシャー情報量は以下で定義される.
$$
J_n(\theta) = E\bigg[ \Big(\frac {\partial} {\partial \theta} \log f(X_1, ..., X_n;\theta) \Big)^2 \bigg]
$$
もしくは、適当な正則条件のもと同値となる
$$
J_n(\theta) = -E\bigg[ \frac {\partial^2} {\partial^2 \theta} \log f(X_1, ..., X_n;\theta) \bigg]
$$
フィッシャー情報量が正であるとき, 不偏推定量に対するクラメール・ラオの不等式は以下で表される.
$$
V[\hat{\theta}] \geq J_n(\theta)^{-1}
$$
これは, 不偏推定量 $\hat{\theta}$ の分散はがフィッシャー情報量の逆数よりも小さくすることはできないことを意味している. なお, 標本が独立であるときは $ J_n(\theta) = n J_1(\theta)$ となりフィッシャー情報量は標本サイズに比例するため, クラメール・ラオの不等式の逆数は標本サイズの逆数に比例する.
十分統計量
以下の式を満たす統計量 $T = T(X)$ を $\theta$ の十分統計量という.
$$
P(X = x | T(X) = t, \theta) = P(X = x | T(X) = t)
$$
つまり, $T(X)$ で条件づけた $X$ の分布がパラメータによらないことが十分統計量の定義である. $T(X)$ が $\theta$ の十分統計量であるとき, 適当な関数 $h,g$ が存在して
$$
f(x,\theta) = h(x)g(T(x), \theta)
$$
という分解が可能であることが示される(フィッシャー・ネイマンの分解定理). つまり, 確率密度関数を $\theta$ に依存しない関数と, 依存する関数の積に分解したときに, 後者が十分統計量のみを含むような分解が存在する.
例4
\begin{align}
f(x;\mu, \sigma^2) &= \prod_{i=1}^{n} \frac {1} {\sqrt{2\pi \sigma^2}} \exp \bigg( - \frac {(x_i - \mu)^2} {2\sigma^2}\bigg) \\
&= (2\pi \sigma^2)^{-n/2} \exp \bigg( - \frac {\sum_{i=1}^{n} x_i^2} {2\sigma^2} \bigg) \exp \bigg( - \frac {n\mu^2} {2\sigma^2} \bigg) \exp \bigg(\frac {n\bar{x}\mu} {\sigma^2} \bigg) \\
&= (2\pi \sigma^2)^{-n/2} \exp \bigg( - \frac {\sum_{i=1}^{n}(x_i - \bar{x})^2} {2\sigma^2} \bigg) \exp \bigg( - \frac {n(\bar{x} - \mu)^2} {2\sigma^2} \bigg)
\end{align}
分散 $\sigma^2$ が既知の時に, 標本平均 $\bar{X}$ は $\mu$ の十分統計量である. これは,
$$
h(x) = (2\pi \sigma^2)^{-n/2} \exp \bigg( - \frac {\sum_{i=1}^{n} x_i^2} {2\sigma^2} \bigg), g(\bar{x}, \mu) = \exp \bigg( - \frac {n\mu^2} {2\sigma^2} \bigg) \exp \bigg(\frac {n\bar{x}\mu} {\sigma^2} \bigg)
$$
とおくことによりフィッシャー・ネイマンの分解定理を用いて確認できる. また, $\mu, \sigma$ がどもに未知のとき, 同じく式に関する分解定理を用いて確認できる. $\bar{x}, \sum_{i=1}^{n} x_i^2$ は $\theta = (\mu, \sigma)$ の十分統計量であることが分かる.
漸近的な性質
標本サイズが十分大きいときの推定の妥当性を評価するための理論を漸近論という.
一致性
推定量 $\hat{\theta}$ が真のパラメータ $\theta$ に確率収束する時, つまり任意の $\epsilon > 0$ に対して,
$$
\lim_{n \to \infty} P(|\hat{\theta} - \theta|) < \epsilon) = 1
$$
が成立する時, 推定量が一致性を持つという.
漸近有効性
一致推定量の分散が漸近的にクラメール・ラオの不等式の下限を達成する時, つまり任意の $\theta$ に対して,
$$
\lim_{n \to \infty} nV[\hat{\theta}] = J_1(\theta)^{-1}
$$
となるとき, 推定量が漸近有効性をもつという.
漸近正規性
最尤推定量 $\hat{\theta}$ から極限値の $\theta$ を引いたうえで $\sqrt{n}$ でスケーリングした
$$
Z = \sqrt{n}(\hat{\theta} - \theta)
$$
は正規分布 $N(0, J_1(\theta)^{-1})$ に分布収束する. これを最尤推定量の漸近正規性という.
リサンプリング法
推定量が不偏でなくバイアスがある場合, 得られている標本を再利用して推定量のバイアスを補正する方法がジャックナイフ法である. 標本 $X_1, ..., X_n$ をもとにした推定量を $\hat{\theta}$ のバイアスを補正する場合を考える.
標本 $X_i$ を除いたサイズ $n-1$ の部分標本から推定量 $\hat{\theta}$ と同様にして計算される推定量を $\hat{\theta}_{(i)}$ と書く. またそのようにして得られた $n$ 個の推定量の平均を
$$
\hat{\theta}_{\cdot} = \frac {1} {n} \sum^{n}_{i=1} \hat{\theta}_{(i)}
$$
とする.
このとき, $\hat{\theta}$ のバイアスは $bias = (n-1)(\hat{\theta}_{(\cdot)} - \hat{\theta})$ で見積もられ, それを補正した
$$
\hat{\theta}_{jack} = \hat{\theta} - bias = \hat{\theta} - (n-1)(\hat{\theta}_{(\cdot)} - \hat{\theta}) = n\hat{\theta} - (n-1) \hat{\theta}_{(\cdot)}
$$
をジャックナイフ推定量という. ジャックナイフ法のように, 得られている標本からの部分標本を持ちリいることにより推定精度を向上させる手法をリサンプリング法と呼ぶ.
例題
問8.2
[1]
$\lambda$ に関する対数尤度関数は,
\begin{align}
l(\lambda) &= \log \prod_{i=1}^{n} e^{-\lambda} \frac {\lambda^{x_i}} {x_i!} \\
&= \sum^{n}_{i=1} \log e^{-\lambda} \frac {\lambda^{x_i}} {x_i!} \\
&= -n\lambda + \sum_{i=1}^{n} x_i \log \lambda - \sum_{i=1}^{n} \log x_i!
\end{align}
これを$\lambda$で微分して,
\begin{align}
\frac {d l(\lambda)} {d \lambda} = -n + \sum_{i=1}^{n} \frac {x_i} {\lambda} = 0
\end{align}
を解いた $\hat{\lambda} = 1/n \sum_{i=1}^{n} x_i$ が最尤唯定量である.
[2]
フィッシャー情報量は,
\begin{align}
J_n(\lambda) &= -E \bigg[ \frac {d^2} {d \lambda^2} (-n + \sum_{i=1}^{n} \frac {x_i} {\lambda}) \bigg] \\
&= -E \bigg[ \frac {d l(\lambda)} {d \lambda} \Big(-n + \sum_{i=1}^{n} \frac {x_i} {\lambda} \Big) \bigg] \\
&= - E \bigg[ \sum_{i=1}^{n} X_i \frac {-1} {\lambda^2} \bigg] \\
&= \frac {1} {\lambda^2} E \bigg[ \sum_{i=1}^{n} X_i \bigg] \\
&= \frac {1} {\lambda^2} \Big( E[X_1] + ... + E[X_n] \Big) \bigg] \\
&= \frac {1} {\lambda^2} n\lambda = \frac {n} {\lambda}
\end{align}
[3]
[1]より, 最尤推定量の分散は
\begin{align}
V \bigg[ \frac {1} {n} \sum_{i=1}^{n} X_i \bigg] &= \frac {1} {n^2} V \bigg[\sum_{i=1}^{n} X_i \bigg] \\
&= \frac {1} {n^2} \Big( V[X_1] + ... + V[X_n] \Big) \\
&= \frac {1} {\lambda^2} n\lambda = \frac {n} {\lambda}
\end{align}
であるが, これはクラメール・ラオの不等式の下限
$$
J_n(\lambda)^{-1} = \frac {\lambda} {n}
$$
と一致する.よって, 標本平均は有効推定量である.
[4]
標本平均に対する中心極限定理より, $\sqrt{n}(\bar{X} - \lambda)$ は $N(0,1)$ に分布収束する. よって, 紺場合の最尤推定量の漸近正規性が示された. また、この場合の漸近分散 $\lambda$ はフィッシャー情報量の逆数 $J_1(\lambda)^{-1} = \lambda$ と一致しており, 漸近有効性も成立している.
問8.3
[1]
各推定量の期待値を計算すると,
\begin{align}
E[T_1] &= E[\bar{X}^2] = E[\bar{X}]^2 + V[X] \\
&= \mu^2 + \frac {1} {n} \sigma^2 \\
E[T_2] &= E\bigg[ \frac {1} {n} \sum_{i=1}^{n} X_i^2 \bigg] \\
&= \frac {1} {n} n E[X_i^2] = E[X_i^2] = \mu^2 + \sigma^2
\end{align}
よって,
\begin{align}
b_1(\mu, \sigma^2) &= E[T_1] - \mu^2 \\
&= \frac {1} {n} \sigma^2 \\
b_2(\mu, \sigma^2) &= E[T_2] - \mu^2 \\
&= \sigma^2
\end{align}
[2]
問題のコインの面積の推定法は, 推定量 $T_2$ を $\pi$ 倍したものに対応する.
[1] の結果より推定量 $T_2$ のバイアスは $\sigma^2$ であり $n$ に依存せず, 標本サイズに関わらないバイアスを持つためこの推定量は問題である.
[3]
[1]より, $b_1(\mu, \sigma^2) = \sigma^2/n $ なので, $\bar{T_1} = \bar{X}^2 - U/n$ である.
5 つの半径の標本平均 $\bar{X}$ は9.04, 不偏分散 $U$ は0.748なので, コインの面積の推定値は
$$
3.14 \times T_1 = 3.14 \times (9.04^2 - 0.748/5) = 256
$$ である.
[4]
$i$ 番目のデータを除いた残りの 4 つで $T_1$ を計算した値を $T_{1(i)}$ とすると, $T_{1(1)} = {(9.0+8.2+10.4+9.2)/4}^2 = 84.64$ となる. 同様にして, $T_{1(2)} = 81.90, T_{1(3)} = 85.56, T_{1(4)} = 75.69, T_{1(5)} = 81.00$ と計算できる. これらの平均を $T_{1(\cdot)} = (T_{1(1)} + T_{1(2)} + T_{1(3)} + T_{1(4)} + T_{1(5)})/5 = 81.76$ としたとき, $T_1$ のジャックナイフ推定量は,
$$
T_{jack} = nT_1 - (n-1)T_{1(\cdot)} = 5 \times 9.04^2 - 4 \times 81.76 = 81.57
$$
よって, コインの面積のバイアス補正推定値は, $3.14 \times T_{jack} = 3.14 \times 81.76 = 256$ となる.
[5]
$\bar{X}$ についての中心極限定理より, $\sqrt{n}(\bar{X} - \mu)$ は $N(0,\mu)$ に分布収束する. さらに, $g(x) = x^2$ という変数変換に対してデルタ法を適用すると, 漸近分散は $\sigma^2 g^{'}(\mu)^2 = 4\mu^2\sigma^2$ となり, $\sqrt{n}(\bar{X} - \mu)$ は$N(0,4\mu^2\sigma^2)$ に分布収束する. よって, 漸近正規性が成立する.