第六章 統計的推定

この章で目指すこと

これまでの章では、確率分布が分かっているときに、平均、分散、確率、標本分布を計算してきました。しかし実際の研究では、分布の形を決める母数は未知です。手元にあるのは有限個のデータだけです。

統計的推定とは、観測データから未知母数について学ぶための方法です。

例えば、薬物投与後の反応量を正規分布で表すとき、母平均 $\mu$ や母分散 $\sigma^2$ は見えません。観測されたマウスや患者の値から、それらを推定します。

手元にある情報	推定したいもの	主な方法	何が学べるか
成功・失敗の回数	反応確率 $p$	最尤法、Beta-Binomialベイズ	真の反応率と不確実性
単位時間のイベント数	発生率 $\lambda$	Poisson最尤、Gamma-Poissonベイズ	有害事象率、コロニー数、発現率
連続測定値	平均 $\mu$、分散 $\sigma^2$	モーメント法、最尤法	平均反応と個体差
複数施設・複数個体の推定値	共通傾向と群ごとの差	階層・経験Bayes	情報共有後の群別推定
尤度曲線の曲がり方	推定精度	Fisher情報量	標準誤差の理論的限界

この章を通して使う小さな例

架空の薬効スクリーニングで、独立な10個体のうち4個体が反応したとします。反応を1、非反応を0とすれば、データの一例は

$$(1,0,1,0,0,1,0,0,1,0), \qquad n=10,\quad S=\sum_{i=1}^{10}x_i=4$$

です。未知なのは、同じ条件での真の反応確率 $p$ です。この1組のデータから、後の節で次を順番に導きます。

見方	この例で得られるもの
十分統計量	10個の並びを成功回数 $S=4$ に圧縮できる
最尤法	$\hat p_{\mathrm{ML}}=S/n=0.4$
Bayes法	$\operatorname{Beta}(2,2)$ 事前分布なら事後分布は $\operatorname{Beta}(6,8)$
Fisher情報量	$p$ の近くで尤度がどれほど鋭く曲がるかを数値化する
Cramér–Rao下限	不偏推定量の分散がどこまで小さくなり得るかを調べる

1. 推定の言葉を整理する

1.1 統計モデル、母数、標本

確率密度または確率関数を $f(x\mid\theta)$ と書きます。$\theta$ は未知母数です。

$$X_1,\ldots,X_n\overset{\mathrm{iid}}{\sim}f(x\mid\theta)$$

例えばBernoulli分布なら $\theta=p$、正規分布なら $\theta=(\mu,\sigma^2)$ です。

• 推定量 $\hat\theta=T(X_1,\ldots,X_n)$：観測前には確率変数
• 推定値 $\hat\theta_{\mathrm{obs}}=T(x_1,\ldots,x_n)$：実データから得た数値
• 母数 $\theta$：固定されているが未知の値

通し例では、観測前の成功回数 $S=\sum_iX_i$ は確率変数なので、

$$\hat p=\frac{S}{10}$$

も確率変数、すなわち推定量です。実際に $S=4$ を観測した後の $0.4$ は推定値です。一方、真の $p$ は観測後も未知の母数です。

記号	観測前	観測後	通し例
$p$	固定だが未知	固定だが未知	真の反応確率
$\hat p=S/10$	標本ごとに変わる確率変数	計算式	推定量
$\hat p_{\mathrm{obs}}$	まだ存在しない	具体的な数値	$4/10=0.4$

考え方	母数 $\theta$	データ $X$	不確実性
頻度論	固定だが未知	繰り返すと変わる	推定量の標本分布
Bayes	事前分布を置く	観測後は固定された情報	事後分布

1.2 点推定と損失

点推定は未知母数を1つの数で表します。ただし、何を「良い推定」とするかには基準が必要です。推定値 $a$ と真値 $\theta$ のずれを損失 $L(\theta,a)$ で表します。

損失	式	Bayes推定量
二乗誤差	$(a-\theta)^2$	事後平均
絶対誤差	$\lvert a-\theta\rvert$	事後中央値
0–1損失の近似	外れたら1	事後最頻値 MAP

毒性の見逃しと偽陽性の損失は同じとは限りません。応用では「何を小さくしたいか」を先に決めます。

2. 十分統計量

2.1 母数の情報を失わない圧縮

標本全体から統計量 $T$ だけを残しても、$\theta$ を推定する情報が失われないなら、$T$ を十分統計量といいます。

Bernoulli標本では、成功の順序ではなく成功回数

$$S=\sum_{i=1}^nX_i$$

が $p$ の情報を持ちます。

例えば、次の2つのデータは並び順が異なります。

$$\begin{aligned} x&=(1,0,1,0,0,1,0,0,1,0),\\ y&=(0,1,0,1,0,0,1,0,0,1). \end{aligned}$$

しかし、どちらも成功4回、失敗6回です。独立性より、各データが現れる確率は

$$p\,(1-p)\,p\,(1-p)\,(1-p)\,p\,(1-p)\,(1-p)\,p\,(1-p) =p^4(1-p)^6$$

となります。$p$ を推定する限り、成功が何番目に現れたかではなく、成功と失敗の個数が重要です。これが「情報を失わない圧縮」の最初のイメージです。

2.2 定義と因子分解定理

$T=t$ を与えた条件付き分布 $\mathcal{L}(X\mid T=t)$ が $\theta$ に依存しないとき、$T$ は十分統計量です。

実際の判定にはNeyman–Fisherの因子分解定理を使います。同時密度または同時確率関数が

$$f_\theta(x_1,\ldots,x_n) =g_\theta\{T(x_1,\ldots,x_n)\}\,h(x_1,\ldots,x_n)$$

と分解できれば、$T$ は十分統計量です。

Bernoulli標本

$$\begin{aligned} L(p) &=\prod_{i=1}^n p^{x_i}(1-p)^{1-x_i}\\ &=p^{\sum_i x_i}(1-p)^{n-\sum_i x_i}. \end{aligned}$$

$p$ を含む部分は $\sum_i x_i$ だけを通してデータに依存します。したがって $T=\sum_iX_i$ は $p$ の十分統計量です。

定義からも確かめられます。$S=4$ と分かった後、成功4個を10か所へ配置する並びは $\binom{10}{4}$ 通りです。どの並びも確率 $p^4(1-p)^6$ を持つので、特定の並び $x$ の条件付き確率は

$$\begin{aligned} P_p(X=x\mid S=4) &=\frac{P_p(X=x)}{P_p(S=4)}\\ &=\frac{p^4(1-p)^6} {\binom{10}{4}p^4(1-p)^6}\\ &=\frac1{\binom{10}{4}}. \end{aligned}$$

最後の式から $p$ が消えました。つまり、$S=4$ を知った後に残る「並び順」は、$p$ について追加情報を持ちません。

支持範囲に母数が入る例

$X_i\sim U(0,\theta)$ なら、

$$L(\theta)=\theta^{-n}I(0<x_{(n)}<\theta), \qquad x_{(n)}=\max_i x_i.$$

$\theta$ を含む部分は最大値を通じてデータを見ているので、$X_{(n)}$ は十分統計量です。

因子分解定理を条件付き分布から導く

離散型で $T(X)=t$ とします。因子分解できるなら、

$$P_\theta(X=x\mid T=t) =\frac{g_\theta(t)h(x)} {\sum_{y:T(y)=t}g_\theta(t)h(y)} =\frac{h(x)}{\sum_{y:T(y)=t}h(y)}.$$

$g_\theta(t)$ が約分され、条件付き分布から $\theta$ が消えます。$T$ を知った後の残りの並び方は、母数について追加情報を持ちません。

2.3 最小十分統計量

十分統計量の中でも、それ以上は本質的に圧縮できないものを最小十分統計量といいます。標本点 $x,y$ に対する尤度比 $f_\theta(x)/f_\theta(y)$ が $\theta$ に依存しないことと、$T(x)=T(y)$ が同値になるかを調べる方法がよく使われます。

Interactive estimation lab

推定量の振る舞いを動かして確かめる

タブを切り替え、標本数や観測値を動かしてください。線や数値はすべてその場で再計算されます。

十分統計量最尤とベイズバイアスとMSE情報量と漸近分布

試行回数 n

8

成功回数 s

3

並び順を捨てても尤度は変わらない

ベルヌーイ標本では、成功の順序ではなく合計 S = ΣXᵢ だけが p の情報を持ちます。

S = s の尤度 最尤推定値

3. 指数型分布族

3.1 一般形

指数型分布族を学ぶ目的は、分布名を新しく覚えることではありません。密度を共通の形へ書き直すと、次の3つをまとめて読めることが利点です。

1. データのどの要約が十分統計量か
2. 尤度を微分したとき何が現れるか
3. 平均・分散が対数分配関数 $A$ の微分からどう得られるか

$$f(x\mid\theta) =h(x)\exp\left\{ \eta(\theta)^\mathsf{T}T(x)-A(\theta) \right\}$$

と書ける分布族を指数型分布族といいます。

• $h(x)$：母数を含まない部分
• $T(x)$：十分統計量の核
• $\eta(\theta)$：自然母数
• $A(\theta)$：対数分配関数

$n$ 個のiid標本では、

$$\prod_{i=1}^nf(x_i\mid\theta) =\left\{\prod_{i=1}^nh(x_i)\right\} \exp\left[ \eta(\theta)^\mathsf{T}\sum_{i=1}^nT(x_i)-nA(\theta) \right].$$

したがって $\sum_iT(X_i)$ が十分統計量になります。

ベクトル式を1母数の式に戻す

$\eta^\mathsf{T}T$ は内積です。1母数なら単に $\eta(\theta)T(x)$、2母数なら $\eta_1T_1+\eta_2T_2$ です。

3.2 Bernoulli、Poisson、正規分布

Bernoulli分布を1行ずつ変形する

$x$ は0か1なので、確率関数は $p^x(1-p)^{1-x}$ と1本に書けます。指数型へ直すため、$a=\exp(\log a)$ を使います。

$$\begin{aligned} f(x\mid p) &=p^x(1-p)^{1-x}\\ &=\exp\{x\log p+(1-x)\log(1-p)\}\\ &=\exp\left[ x\{\log p-\log(1-p)\}+\log(1-p) \right]\\ &=\exp\left[ x\log\frac{p}{1-p}+\log(1-p) \right]. \end{aligned}$$

したがって自然母数は $\eta=\log\{p/(1-p)\}$、統計量は $T(x)=x$ です。自然母数は確率 $p$ をlogit変換した量で、$p=e^\eta/(1+e^\eta)$ と元へ戻せます。

Poisson分布

Poisson分布は

$$f(x\mid\lambda)=\frac1{x!}\exp\{x\log\lambda-\lambda\}$$

なので、自然母数は $\eta=\log\lambda$ です。

正規分布で平方を展開する

正規分布では、指数部の平方を

$$(x-\mu)^2=x^2-2\mu x+\mu^2$$

と展開します。これを密度へ戻すと、

$$\begin{aligned} f(x\mid\mu,\sigma^2) &=\frac1{\sqrt{2\pi\sigma^2}} \exp\left\{-\frac{(x-\mu)^2}{2\sigma^2}\right\}\\ &=\exp\left[ \frac{\mu}{\sigma^2}x-\frac{x^2}{2\sigma^2} -\frac{\mu^2}{2\sigma^2}-\frac12\log(2\pi\sigma^2) \right]. \end{aligned}$$

したがって、

$$T(x)=\begin{pmatrix}x\\x^2\end{pmatrix}, \qquad \eta=\begin{pmatrix}\mu/\sigma^2\\-1/(2\sigma^2)\end{pmatrix}.$$

標本では $\left(\sum_iX_i,\sum_iX_i^2\right)$ が十分統計量です。

分布	自然母数 $\eta$	$T(x)$	期待値母数
Bernoulli	$\log\{p/(1-p)\}$	$x$	$E[X]=p$
Poisson	$\log\lambda$	$x$	$E[X]=\lambda$
正規、$\sigma^2$ 既知	$\mu/\sigma^2$	$x$	$E[X]=\mu$
正規、両方未知	$(\mu/\sigma^2,-1/(2\sigma^2))^\mathsf T$	$(x,x^2)^\mathsf T$	$(E[X],E[X^2])^\mathsf T$

この2次元ベクトルで実際に何を計算しているか

内積は $(\mu/\sigma^2)x+\{-1/(2\sigma^2)\}x^2$ です。観測値の和と二乗和を保存すれば、$\mu$ と $\sigma^2$ の尤度を再構成できます。$\sigma^2$ が既知なら必要なのは $\sum X_i$ だけです。

3.3 自然母数と期待値母数

まず2つの母数の役割を分ける

自然母数表示

$$f(x\mid\eta) =h(x)\exp\{\eta^\mathsf{T}T(x)-A(\eta)\}$$

では、同じ分布を2つの方向から表せます。

表し方	定義	何を表すか	便利な場面
自然母数 $\eta$	指数部で $T(x)$ に掛かる係数	尤度がデータをどの方向へ重視するか	尤度の微分、最適化、GLM
期待値母数 $m$	$m=E_\eta[T(X)]$	十分統計量が平均的にどの値を取るか	解釈、モーメント、最尤方程式

重要なのは、期待値母数が常に $E[X]$ とは限らないことです。$T(X)=X$ なら $E[X]$ ですが、$T(X)=(X,X^2)^\mathsf T$ なら

$$m= \begin{pmatrix} E[X]\\ E[X^2] \end{pmatrix}$$

です。

自然母数は、元の母数を指数部の係数として最も整理しやすい形へ変換したものです。例えばBernoulli分布では元の母数が $p$、自然母数が

$$\eta=\log\frac{p}{1-p}$$

です。$p$ と $\eta$ は同じ分布を指していますが、目盛りが違います。

$p$	自然母数 $\eta=\log\{p/(1-p)\}$	読み方
0.2	$-1.386$	成功より失敗が起こりやすい
0.5	$0$	成功と失敗が同程度
0.8	$1.386$	成功の方が起こりやすい

$A(\eta)$ は何をしているか

$A(\eta)$ は対数分配関数です。確率の総和または密度の積分を1にするための正規化定数を、対数で表しています。

1母数の連続型なら、

$$\begin{aligned} 1 &=\int f(x\mid\eta)\,dx\\ &=\int h(x)\exp\{\eta T(x)-A(\eta)\}\,dx\\ &=e^{-A(\eta)} \int h(x)e^{\eta T(x)}\,dx. \end{aligned}$$

したがって、

$$e^{A(\eta)} =\int h(x)e^{\eta T(x)}\,dx$$

であり、

$$A(\eta) =\log\int h(x)e^{\eta T(x)}\,dx$$

です。離散型では積分を総和へ置き換えます。つまり、$\eta$ を決めると $A(\eta)$ も自動的に決まり、密度全体が1になるよう調整されます。

1階微分から期待値が出る導出

正規化条件

$$\int f(x\mid\eta)\,dx=1$$

を $\eta$ で微分します。右辺の微分は0です。一方、

$$\frac{\partial}{\partial\eta}f(x\mid\eta) =\{T(x)-A'(\eta)\}f(x\mid\eta)$$

なので、

$$\begin{aligned} 0 &=\frac{\partial}{\partial\eta} \int f(x\mid\eta)\,dx\\ &=\int\{T(x)-A'(\eta)\}f(x\mid\eta)\,dx\\ &=\int T(x)f(x\mid\eta)\,dx -A'(\eta)\int f(x\mid\eta)\,dx\\ &=E_\eta[T(X)]-A'(\eta). \end{aligned}$$

よって、

$$A'(\eta)=E_\eta[T(X)]$$

です。対数分配関数の傾きが期待値母数になります。

2階微分から分散が出る導出

もう一度微分すると、

$$\begin{aligned} A''(\eta) &=\frac{\partial}{\partial\eta}E_\eta[T(X)]\\ &=\int T(x)\frac{\partial f(x\mid\eta)} {\partial\eta}\,dx\\ &=\int T(x)\{T(x)-A'(\eta)\} f(x\mid\eta)\,dx\\ &=E_\eta[T(X)^2] -E_\eta[T(X)]^2\\ &=\operatorname{Var}_\eta\{T(X)\}. \end{aligned}$$

したがって、

$$A'(\eta)=E_\eta[T(X)], \qquad A''(\eta)=\operatorname{Var}_\eta\{T(X)\}.$$

分散は0以上なので $A''(\eta)\ge0$ です。したがって $A$ は凸関数です。分散が正なら $A'$ は単調増加し、自然母数 $\eta$ と期待値母数 $m=A'(\eta)$ を1対1に行き来できます。

Bernoulli分布

Bernoulli分布では、

$$\eta=\log\frac{p}{1-p}, \qquad p=\frac{e^\eta}{1+e^\eta}$$

です。自然母数表示は

$$f(x\mid\eta) =\exp\left\{x\eta-\log(1+e^\eta)\right\}$$

なので、

$$A(\eta)=\log(1+e^\eta).$$

1階微分は、

$$\begin{aligned} A'(\eta) &=\frac{e^\eta}{1+e^\eta}\\ &=p\\ &=E[X]. \end{aligned}$$

2階微分は、

$$\begin{aligned} A''(\eta) &=\frac{e^\eta(1+e^\eta)-e^{2\eta}} {(1+e^\eta)^2}\\ &=\frac{e^\eta}{(1+e^\eta)^2}\\ &=p(1-p)\\ &=\operatorname{Var}(X). \end{aligned}$$

ここでは $T(X)=X$ なので、期待値母数は $m=p$ です。

Poisson分布

Poisson分布では、

$$\eta=\log\lambda, \qquad \lambda=e^\eta$$

であり、

$$f(x\mid\eta) =\frac1{x!}\exp\{x\eta-e^\eta\}.$$

したがって、

$$A(\eta)=e^\eta, \qquad A'(\eta)=e^\eta=\lambda, \qquad A''(\eta)=e^\eta=\lambda.$$

Poisson分布で平均と分散がともに $\lambda$ になることが、$A$ の1階・2階微分から同時に出ます。

正規分布：分散 $\sigma^2$ が既知の場合

$\sigma^2$ を固定すると、

$$\eta=\frac{\mu}{\sigma^2}, \qquad A(\eta)=\frac{\sigma^2\eta^2}{2}$$

と書けます。したがって、

$$A'(\eta)=\sigma^2\eta=\mu, \qquad A''(\eta)=\sigma^2.$$

ここでも $T(X)=X$ なので、期待値母数は $E[X]=\mu$、その変化率は分散 $\sigma^2$ です。

正規分布：平均と分散がともに未知の場合

この場合は、

$$T(X)= \begin{pmatrix} X\\X^2 \end{pmatrix}, \qquad \eta= \begin{pmatrix} \eta_1\\\eta_2 \end{pmatrix} = \begin{pmatrix} \mu/\sigma^2\\-1/(2\sigma^2) \end{pmatrix}.$$

したがって期待値母数は、

$$m= \begin{pmatrix} E[X]\\E[X^2] \end{pmatrix} = \begin{pmatrix} \mu\\\mu^2+\sigma^2 \end{pmatrix}.$$

第2成分は分散ではなく2次モーメントです。期待値母数から元の母数へ戻すときは、

$$\mu=m_1, \qquad \sigma^2=m_2-m_1^2$$

と計算します。これは $\operatorname{Var}(X)=E[X^2]-E[X]^2$ そのものです。

自然母数と期待値母数では、許される範囲も違って見えます。

分布	自然母数の範囲	期待値母数の範囲
Bernoulli	$\eta\in\mathbb R$	$0<m<1$
Poisson	$\eta\in\mathbb R$	$m>0$
正規、両方未知	$\eta_1\in\mathbb R,\ \eta_2<0$	$m_1\in\mathbb R,\ m_2>m_1^2$

正規分布で $m_2>m_1^2$ が必要なのは、$\sigma^2=m_2-m_1^2>0$ でなければならないためです。

正規分布の $A(\eta)$ を偏微分して確かめる

自然母数空間は $\eta_2<0$ で、$A(\eta)=-\eta_1^2/(4\eta_2)+\frac12\log\{\pi/(-\eta_2)\}$ と書けます。偏微分すると $\partial A/\partial\eta_1=-\eta_1/(2\eta_2)=\mu$、$\partial A/\partial\eta_2=\eta_1^2/(4\eta_2^2)-1/(2\eta_2)=\mu^2+\sigma^2$ となります。2変数でも「$A$ の各方向の傾き＝対応する十分統計量の期待値」です。

3.4 最尤法では「標本の値＝期待値母数」になる

$n$ 個のiid標本の対数尤度を、自然母数の関数として書くと、

$$\ell(\eta) =\eta^\mathsf T\sum_{i=1}^nT(x_i) -nA(\eta) +\sum_{i=1}^n\log h(x_i).$$

$\eta$ で微分すると、

$$\nabla_\eta\ell(\eta) =\sum_{i=1}^nT(x_i) -n\nabla A(\eta).$$

内点の最尤推定値ではスコアを0と置くので、

$$\frac1n\sum_{i=1}^nT(x_i) =\nabla A(\hat\eta) =E_{\hat\eta}[T(X)].$$

つまり最尤法は、観測された十分統計量の平均とモデルが予測する十分統計量の平均が一致するように母数を選びます。

分布	左辺：標本から得る量	右辺：モデルの期待値	得られる最尤推定値
Bernoulli	$\bar X$	$p$	$\hat p=\bar X$
Poisson	$\bar X$	$\lambda$	$\hat\lambda=\bar X$
正規、両方未知	$(\bar X,\overline{X^2})$	$(\mu,\mu^2+\sigma^2)$	$\hat\mu=\bar X$、$\hat\sigma^2=\overline{X^2}-\bar X^2$

自然母数表示を使うと、この3つの最尤計算が同じ1本の式で説明できます。

多次元では勾配とHessianになる

$$\nabla A(\eta)=E_\eta[T(X)], \qquad \nabla^2A(\eta)=\operatorname{Cov}_\eta\{T(X)\}.$$

Hessianは共分散行列なので半正定値であり、$A(\eta)$ は凸関数です。

1変数なら何に当たるか

勾配は1階微分 $A'(\eta)$、Hessianは2階微分 $A''(\eta)$、共分散行列は分散 $\operatorname{Var}(T)$ に戻ります。

発展：自然母数と期待値母数は凸双対の座標

$A$ が狭義凸なら $m=\nabla A(\eta)$ は1対1で、逆向きは凸共役 $A^*(m)=\sup_\eta\{\eta^\mathsf Tm-A(\eta)\}$ から $\eta=\nabla A^*(m)$ と表せます。自然母数は尤度計算に適した座標、期待値母数は分布の平均的な振る舞いを読む座標です。これは情報幾何や最大エントロピー法へつながりますが、統計検定1級では、まず $\nabla A=E[T]$ と $\nabla^2A=\operatorname{Cov}(T)$ を使えることが中心です。

4. 推定量を比べる基準

4.1 バイアス、不偏性、平均二乗誤差

推定量 $\hat\theta$ のバイアスは

$$\operatorname{Bias}_\theta(\hat\theta) =E_\theta[\hat\theta]-\theta$$

です。すべての $\theta$ で0なら不偏推定量です。

平均二乗誤差は

$$\operatorname{MSE}_\theta(\hat\theta) =E_\theta[(\hat\theta-\theta)^2]$$

であり、

$$\operatorname{MSE}(\hat\theta) =\operatorname{Var}(\hat\theta) +\operatorname{Bias}(\hat\theta)^2$$

と分解できます。

MSE分解の途中計算

$$\hat\theta-\theta =(\hat\theta-E\hat\theta)+(E\hat\theta-\theta)$$

を二乗して期待値を取ると、

$$\begin{aligned} E[(\hat\theta-\theta)^2] &=E[(\hat\theta-E\hat\theta)^2]\\ &\quad+2(E\hat\theta-\theta)E[\hat\theta-E\hat\theta]\\ &\quad+(E\hat\theta-\theta)^2. \end{aligned}$$

$E[\hat\theta-E\hat\theta]=0$ なので中央の項が消えます。

不偏であるだけでは最良とは限りません。機械学習の正則化やBayesの縮小推定は、少しバイアスを許して分散を下げ、MSEや予測誤差を改善します。

数値で比べる

真値を $\theta=10$ とし、2つの推定量の性質が次の通りだったとします。

推定量	$E[\hat\theta]$	Bias	分散	MSE
A	10	0	4	$4+0^2=4$
B	11	1	1	$1+1^2=2$

Aは不偏ですが、ばらつきが大きいためMSEは4です。Bは1だけ上に偏るものの、ばらつきが小さくMSEは2です。「不偏だから常に良い」ではなく、バイアスと分散の合計で目的に合うかを見ます。

4.2 一致性

標本数を増やすと推定量が真値へ確率収束するとき、

$$\hat\theta_n\xrightarrow{p}\theta,$$

$\hat\theta_n$ は一致推定量です。

不偏性は有限標本での平均的なずれ、一致性は $n\to\infty$ での性質です。どちらか一方から他方は自動的には導けません。

4.3 BLUE

BLUEは、同じものを複数回・複数施設で測ったとき、精度の良い測定を少し強く信頼して1つの値へまとめる方法です。

名前は長いですが、次の4語に分ければ十分です。

文字	意味	この章での読み方
B: Best	いちばん良い	後で決める候補の中で、ばらつきが最小
L: Linear	線形	重み付き平均 $a_1Y_1+\cdots+a_kY_k$ の形
U: Unbiased	不偏	平均的には真値を過大・過小評価しない
E: Estimator	推定量	データから真値を推定する計算式

まずは「なぜ単純平均では足りないのか」

同じ標準試料の真の濃度を $\mu$ とし、2施設が測定した値を

$$Y_1=9.8,\qquad Y_2=10.4$$

とします。どちらの施設も、長く繰り返せば平均的には真値を返す、すなわち

$$E[Y_1]=E[Y_2]=\mu$$

と仮定します。ただし測定の安定さは異なり、

$$\operatorname{Var}(Y_1)=1,\qquad \operatorname{Var}(Y_2)=4$$

とします。施設2は標準偏差が2、施設1は標準偏差が1なので、施設2の方が測定値が散らばりやすい設定です。

単純平均なら

$$\frac{Y_1+Y_2}{2} =\frac{9.8+10.4}{2} =10.1$$

です。しかし、精度が違う2つの測定を50%ずつ信頼してよいとは限りません。BLUEは「より安定な施設1を、どの程度強く重み付けすれば最も安定な推定値になるか」を決めます。

重み付き平均を作る

2つの測定を、

$$\hat\mu(a)=aY_1+(1-a)Y_2$$

と混ぜます。$a=0.5$ なら単純平均、$a=1$ なら施設1だけ、$a=0$ なら施設2だけを使うことです。

重みを $a$ と $1-a$ にしておくと、必ず和が1になります。そのため、

$$\begin{aligned} E[\hat\mu(a)] &=aE[Y_1]+(1-a)E[Y_2]\\ &=a\mu+(1-a)\mu\\ &=\mu. \end{aligned}$$

つまり、どの $a$ を選んでも不偏です。ここまでが U: Unbiased です。

次に、$Y_1$ と $Y_2$ が独立なら、重み付き平均の分散は

$$\begin{aligned} \operatorname{Var}\{\hat\mu(a)\} &=a^2\operatorname{Var}(Y_1) +(1-a)^2\operatorname{Var}(Y_2)\\ &=a^2+4(1-a)^2\\ &=5a^2-8a+4. \end{aligned}$$

となります。重みに2乗が付くのは、定数倍の分散が $\operatorname{Var}(aY)=a^2\operatorname{Var}(Y)$ となるためです。

どの重みが最も安定か

分散 $V(a)=5a^2-8a+4$ を最小にします。

$$\begin{aligned} V'(a)&=10a-8,\\ 10a-8&=0,\\ a&=0.8. \end{aligned}$$

したがって、施設1に80%、施設2に20%の重みを置くと最も安定です。

施設1への重み $a$	何をしているか	分散 $V(a)$
0	施設2だけを使う	4
0.5	単純平均	1.25
0.8	BLUE	0.8
1	施設1だけを使う	1

BLUEの分散0.8は、精度の良い施設1だけを使った場合の1よりも小さくなります。精度の低い施設2でも、独立な測定として少しだけ情報を足せるからです。

実際の統合値は、

$$\begin{aligned} \hat\mu_{\mathrm{BLUE}} &=0.8\times9.8+0.2\times10.4\\ &=9.92. \end{aligned}$$

です。

独立な推定量 $Y_i$ が

$$E[Y_i]=\mu,\qquad \operatorname{Var}(Y_i)=\sigma_i^2$$

を満たすとします。一般に $\hat\mu=\sum_i a_iY_i$ が不偏となる条件は $\sum_i a_i=1$ です。分散

$$\operatorname{Var}(\hat\mu)=\sum_i a_i^2\sigma_i^2$$

をこの条件下で最小化すると、

$$a_i=\frac{\sigma_i^{-2}}{\sum_j\sigma_j^{-2}}.$$

これは逆分散重みです。分散が小さい、つまり測定が安定なものほど、大きな重みを受け取ります。先ほどの例なら、

$$\frac1{\sigma_1^2}:\frac1{\sigma_2^2} =\frac11:\frac14 =4:1$$

なので、正規化した重みは $4/(4+1)=0.8$ と $1/(4+1)=0.2$ です。2つなら、

$$\hat\mu =\frac{Y_1/\sigma_1^2+Y_2/\sigma_2^2} {1/\sigma_1^2+1/\sigma_2^2}.$$

状況	まず使う考え方
全て同じ精度	普通の平均
精度だけが異なる	逆分散重みBLUE
施設・個体ごとに真の差もある	ランダム効果・階層モデル
測定値どうしが相関する	共分散行列を含む一般化BLUE

施設ごとに精度が異なる濃度測定や固定効果メタ解析に現れます。医薬研究では、同じ標準物質の施設間測定を統合する場面が分かりやすい例です。一方、施設間の真の差まである場合は、単純なBLUEではなくランダム効果や階層モデルを検討します。

Gauss–Markov定理の行列表現

$$Y=X\beta+\varepsilon,\qquad E[\varepsilon]=0,\qquad \operatorname{Cov}(\varepsilon)=\sigma^2I$$

では、

$$\hat\beta=(X^\mathsf{T}X)^{-1}X^\mathsf{T}Y$$

が線形不偏推定量の中で共分散行列最小です。

単回帰なら何をしている式か

$X$ は切片の列と説明変数の列を並べた表、$\beta$ は切片と傾きです。説明変数が1つなら、傾きは $\sum(x_i-\bar x)(y_i-\bar y)/\sum(x_i-\bar x)^2$ に戻ります。

5. モーメント法

5.1 基本手順

母集団の $k$ 次モーメント

$$\mu_k'=E_\theta[X^k]$$

と標本モーメント

$$m_k'=\frac1n\sum_{i=1}^nX_i^k$$

を等しいと置き、未知母数について解きます。未知母数が $r$ 個なら、通常は $r$ 本の方程式を使います。

Poisson分布では $E[X]=\lambda$ なので、

$$\hat\lambda_{\mathrm{MM}}=\bar X.$$

例えば5区画のコロニー数が $(2,0,3,1,4)$ なら、

$$\bar x=\frac{2+0+3+1+4}{5}=2$$

なので、モーメント推定値は $\hat\lambda_{\mathrm{MM}}=2$ 個/区画です。ここで行ったことは「理論上の平均 $\lambda$」を「手元の平均2」で置き換えただけです。

5.2 Gamma分布の2母数推定

形状 $\alpha$、尺度 $\beta$ のGamma分布では、

$$E[X]=\alpha\beta,\qquad \operatorname{Var}(X)=\alpha\beta^2.$$

分母を $n$ とした標本2次中心モーメント

$$m_2=\frac1n\sum_{i=1}^n(X_i-\bar X)^2$$

を使い、

$$\bar X=\alpha\beta,\qquad m_2=\alpha\beta^2$$

と置きます。第2式を第1式で割ると、

$$\hat\beta_{\mathrm{MM}}=\frac{m_2}{\bar X}, \qquad \hat\alpha_{\mathrm{MM}}=\frac{\bar X^2}{m_2}.$$

例えば $\bar x=6$、$m_2=18$ なら、まず

$$\hat\beta=\frac{18}{6}=3$$

を求め、次に

$$\hat\alpha=\frac{6}{3}=2 \qquad \left(\text{または }\frac{6^2}{18}=2\right)$$

を得ます。2本の式を同時に眺めるより、割り算で一方の母数を消すと計算が軽くなります。

モーメント法は計算しやすい一方、推定値が母数空間の外へ出たり、最尤推定より効率が低くなったりします。最尤計算の初期値としても使われます。

6. 最尤法

6.1 尤度の意味

観測値 $x=(x_1,\ldots,x_n)$ を固定し、母数の関数として

$$L(\theta\mid x)=\prod_{i=1}^nf(x_i\mid\theta)$$

を見ます。確率は「母数を固定してデータを動かす」、尤度は「データを固定して母数を比べる」ものです。

$$\hat\theta_{\mathrm{ML}} =\arg\max_{\theta\in\Theta}L(\theta\mid x).$$

積を和に変えるため、通常は対数尤度 $\ell(\theta)=\log L(\theta)$ を最大化します。

6.2 答案テンプレート

1. 支持範囲を含めて尤度を書く
2. 対数尤度を取る
3. 内点なら微分してスコア方程式を解く
4. 2階微分、単調性、端点から最大であることを確認する
5. 母数空間に入っているか確認する

6.3 Bernoulli分布

$S=\sum_iX_i$ とすると、

$$\ell(p)=S\log p+(n-S)\log(1-p).$$ $$\ell'(p)=\frac{S}{p}-\frac{n-S}{1-p}=0$$

より、

$$S(1-p)-(n-S)p=0, \qquad \hat p_{\mathrm{ML}}=\frac{S}{n}.$$

さらに

$$\ell''(p)=-\frac{S}{p^2}-\frac{n-S}{(1-p)^2}<0$$

なので内点では最大です。$S=0,n$ なら境界も確認します。

通し例：10個中4個なら

通し例では $n=10$、$S=4$ なので、

$$L(p)=p^4(1-p)^6, \qquad \ell(p)=4\log p+6\log(1-p).$$

微分すると、

$$\ell'(p)=\frac4p-\frac6{1-p}.$$

分母を払って0と置けば、

$$\begin{aligned} \frac4p-\frac6{1-p}&=0\\ 4(1-p)-6p&=0\\ 4-10p&=0\\ \hat p_{\mathrm{ML}}&=0.4. \end{aligned}$$

「成功回数/総数」が突然現れたのではなく、成功側の傾き $4/p$ と失敗側の傾き $6/(1-p)$ が釣り合う点を解いた結果です。

6.4 Poisson分布

$$\ell(\lambda) =\sum_i\{x_i\log\lambda-\lambda-\log(x_i!)\}.$$ $$\ell'(\lambda)=\frac{\sum_i x_i}{\lambda}-n=0$$

から、

$$\hat\lambda_{\mathrm{ML}}=\bar X.$$

先ほどの $(2,0,3,1,4)$ では $\sum_i x_i=10$、$n=5$ なので、

$$\frac{10}{\lambda}-5=0 \ \Longrightarrow\ 10=5\lambda \ \Longrightarrow\ \hat\lambda_{\mathrm{ML}}=2.$$

この例ではモーメント推定量と一致しますが、一般には一致しません。

6.5 正規分布

$$\ell(\mu,\sigma^2) =-\frac n2\log(2\pi)-\frac n2\log\sigma^2 -\frac{1}{2\sigma^2}\sum_i(x_i-\mu)^2.$$

偏微分して、

$$\frac{\partial\ell}{\partial\mu} =\frac{1}{\sigma^2}\sum_i(x_i-\mu)=0$$

より $\hat\mu_{\mathrm{ML}}=\bar X$ です。さらに、

$$\begin{aligned} \frac{\partial\ell}{\partial\sigma^2} &=-\frac{n}{2\sigma^2} +\frac{1}{2\sigma^4}\sum_i(x_i-\mu)^2=0. \end{aligned}$$

両辺を $2\sigma^4$ 倍すると、

$$-n\sigma^2+\sum_i(x_i-\mu)^2=0$$

なので、

$$\sigma^2=\frac1n\sum_i(x_i-\mu)^2.$$

ここへ先に求めた $\hat\mu=\bar x$ を代入して、

$$\hat\sigma^2_{\mathrm{ML}} =\frac1n\sum_i(X_i-\bar X)^2$$

を得ます。

分母が $n-1$ の不偏標本分散ではありません。最尤推定量は必ずしも不偏ではありません。

6.6 一様分布の端点

$$L(\theta)=\theta^{-n}I(\theta\ge x_{(n)}).$$

$\theta\ge x_{(n)}$ では単調減少なので、許される最小値

$$\hat\theta_{\mathrm{ML}}=X_{(n)}$$

が最尤推定量です。微分して0と置く方法では解けません。

6.7 最尤推定量の不変性

$\hat\theta_{\mathrm{ML}}$ が $\theta$ の最尤推定量なら、$g(\theta)$ の最尤推定量は $g(\hat\theta_{\mathrm{ML}})$ です。

例えば $\mathrm{EC}_{50}$ の最尤推定値が得られたなら、$\mathrm{pEC}_{50}=-\log_{10}\mathrm{EC}_{50}$ の最尤推定値はその変換で得られます。ただし標準誤差はデルタ法や尤度区間で評価します。

7. Bayes法

7.1 Bayes更新

$$\pi(\theta\mid x) =\frac{L(\theta\mid x)\pi(\theta)} {\int L(\theta\mid x)\pi(\theta)\,d\theta} \propto L(\theta\mid x)\pi(\theta).$$

• $\pi(\theta)$：事前分布
• $L(\theta\mid x)$：尤度
• $\pi(\theta\mid x)$：事後分布
• 分母：周辺尤度

7.2 Beta–Binomial共役

$$X\mid p\sim\operatorname{Binomial}(n,p), \qquad p\sim\operatorname{Beta}(\alpha,\beta).$$

尤度と事前分布を掛けると、

$$\begin{aligned} \pi(p\mid x) &\propto \underbrace{p^x(1-p)^{n-x}}_{\text{尤度}} \underbrace{p^{\alpha-1}(1-p)^{\beta-1}}_{\text{事前分布}}\\ &=p^{x+\alpha-1}(1-p)^{n-x+\beta-1}\\ &\Rightarrow p\mid x\sim\operatorname{Beta}(\alpha+x,\beta+n-x). \end{aligned}$$

Beta分布の密度は $p^{a-1}(1-p)^{b-1}$ に比例します。したがって、掛け算した後の指数を「$a-1$」「$b-1$」と見比べれば、事後分布の2母数を読めます。

二乗誤差損失のBayes推定量は事後平均なので、

$$\hat p_B=\frac{\alpha+x}{\alpha+\beta+n}.$$

さらに、

$$\hat p_B =\frac{n}{n+\alpha+\beta}\frac{x}{n} +\frac{\alpha+\beta}{n+\alpha+\beta}\frac{\alpha}{\alpha+\beta}.$$

データの割合と事前平均の加重平均になっています。

通し例：事前分布 $\operatorname{Beta}(2,2)$ を使う

10個中4個が反応したので、成功側と失敗側をそれぞれ更新すると、

$$\begin{aligned} p\mid x &\sim\operatorname{Beta}(2+4,\ 2+10-4)\\ &=\operatorname{Beta}(6,8). \end{aligned}$$

事後平均は、

$$E[p\mid x]=\frac6{6+8}=\frac6{14}\approx0.429.$$

加重平均として書けば、

$$\frac{10}{14}\times\underbrace{\frac4{10}}_{\text{データ}} +\frac4{14}\times\underbrace{\frac24}_{\text{事前平均}} =0.429.$$

最尤推定値0.4より少し0.5側へ動いたのは、事前分布が「成功2、失敗2に相当する情報」を加えたためです。また、次の1個体が反応する事後予測確率も

$$P(\tilde X=1\mid x)=E[p\mid x]=\frac6{14}$$

です。

7.3 Gamma–Poisson共役

率母数表示

$$\lambda\sim\operatorname{Gamma}(\alpha,\beta), \qquad \pi(\lambda)\propto\lambda^{\alpha-1}e^{-\beta\lambda}$$

を使います。$X_i\mid\lambda\sim\operatorname{Poisson}(\lambda)$ なら、

$$L(\lambda)\propto\lambda^{\sum_i x_i}e^{-n\lambda},$$

したがって、

$$\lambda\mid x \sim\operatorname{Gamma} \left(\alpha+\sum_i x_i,\ \beta+n\right).$$ $$E[\lambda\mid x] =\frac{\alpha+\sum_i x_i}{\beta+n}, \qquad \operatorname{Var}(\lambda\mid x) =\frac{\alpha+\sum_i x_i}{(\beta+n)^2}.$$

例えば事前分布を率表示の $\operatorname{Gamma}(2,1)$、5区画の合計カウントを8とすると、

$$\lambda\mid x\sim\operatorname{Gamma}(2+8,1+5) =\operatorname{Gamma}(10,6)$$

であり、事後平均は $10/6\approx1.67$ 個/区画です。「形状にはイベント数、率には観測機会の量を足す」と読むと整理しやすくなります。

7.4 事後予測と区間の解釈

新しい観測 $\tilde X$ の予測には、

$$p(\tilde x\mid x) =\int p(\tilde x\mid\theta)\pi(\theta\mid x)\,d\theta$$

を使います。母数の不確実性も積分するため、プラグイン予測より裾が広くなることがあります。

Bayes信用区間は $P(\theta\in C\mid x)=0.95$ と解釈できます。頻度論の95%信頼区間は、同じ手続きを繰り返したとき95%が真値を含む区間構成法です。同じ数値でも意味は異なります。

8. スコア関数とFisher情報量

8.1 スコア関数

対数尤度の1階微分

$$U_n(\theta)=\frac{\partial}{\partial\theta}\ell_n(\theta)$$

をスコア関数といいます。尤度曲線の傾きです。

• $U_n(\theta)>0$：$\theta$ を増やすと尤度が増える
• $U_n(\theta)<0$：$\theta$ を増やすと尤度が減る
• $U_n(\hat\theta)=0$：内点の最尤推定候補

通し例のBernoulli対数尤度では、

$$U_{10}(p)=\frac4p-\frac6{1-p}.$$

例えば、

$$U_{10}(0.2)=20-7.5=12.5>0, \qquad U_{10}(0.7)\approx5.71-20<0.$$

$p=0.2$ では右へ進むと尤度が増え、$p=0.7$ では左へ戻ると尤度が増えます。その間の $p=0.4$ でスコアが0になります。スコアは「推定値そのもの」ではなく、現在の母数をどちらへ動かすべきかを示す傾きです。

8.2 スコアの期待値は0

支持範囲が $\theta$ に依存せず、積分と微分を交換できるとします。1観測のスコアは

$$U_1(\theta)=\frac{\partial}{\partial\theta}\log f(X\mid\theta) =\frac{\partial f(X\mid\theta)/\partial\theta}{f(X\mid\theta)}.$$

したがって、

$$\begin{aligned} E_\theta[U_1(\theta)] &=\int \frac{\partial f(x\mid\theta)/\partial\theta}{f(x\mid\theta)} f(x\mid\theta)\,dx\\ &=\int\frac{\partial}{\partial\theta}f(x\mid\theta)\,dx\\ &=\frac{\partial}{\partial\theta}\int f(x\mid\theta)\,dx\\ &=\frac{\partial}{\partial\theta}1=0. \end{aligned}$$

一様分布 $U(0,\theta)$ のように支持範囲が母数へ依存する場合、この正則条件は破れます。

8.3 Fisher情報量

1観測のFisher情報量を

$$I_1(\theta)=E_\theta[U_1(\theta)^2]$$

で定義します。正則条件のもとでは、

$$I_1(\theta) =-E_\theta\left[ \frac{\partial^2}{\partial\theta^2}\log f(X\mid\theta) \right].$$

2つの式が等しい理由

$E[U]=0$ を $\theta$ で微分すると、密度自体も $\theta$ に依存するため、

$$0=\frac{\partial}{\partial\theta}E[U] =E\left[\frac{\partial U}{\partial\theta}\right]+E[U^2].$$

よって $E[U^2]=-E[\partial U/\partial\theta]$ です。

iid標本ではスコアが和になるので、

$$I_n(\theta)=nI_1(\theta).$$

独立な観測を増やすと情報量は足し算されます。

8.4 代表分布の情報量

Bernoulli分布

$$\begin{aligned} \ell_1(p)&=X\log p+(1-X)\log(1-p),\\ U_1(p)&=\frac{\partial\ell_1}{\partial p}\\ &=\frac{X}{p}-\frac{1-X}{1-p}\\ &=\frac{X-p}{p(1-p)}. \end{aligned}$$

$\operatorname{Var}(X)=p(1-p)$ より、

$$I_1(p) =\frac{\operatorname{Var}(X)}{p^2(1-p)^2} =\frac{1}{p(1-p)}.$$

通し例の推定値 $p=0.4$ を代入すると、

$$I_{10}(0.4) =10\times\frac1{0.4(1-0.4)} \approx41.67.$$

情報量の逆数は $1/41.67\approx0.024$、その平方根は約0.155です。これは標本比率の標準誤差

$$\sqrt{\frac{p(1-p)}n} =\sqrt{\frac{0.4\times0.6}{10}} \approx0.155$$

と一致します。

Poisson分布

ここでは $X\sim\operatorname{Poisson}(\lambda)$、すなわち

$$P(X=x\mid\lambda)=\frac{e^{-\lambda}\lambda^x}{x!}, \qquad x=0,1,2,\ldots$$

とします。$\lambda$ は「一定の観察単位で平均して何件起こるか」を表す母数です。例えば、一定面積の培養皿で生じるコロニー数、一定時間に記録されるカルシウムスパイク数、一定観察時間に起こる有害事象数を表せます。

1観測から、対数尤度をつくる

実現値を小文字 $x$ と書き、$x!$ も含めて $\lambda$ の関数として眺めます。これが尤度です。

$$L_1(\lambda;x)=\frac{e^{-\lambda}\lambda^x}{x!}.$$

対数を取ると、積が和に、べきが係数になります。

$$\begin{aligned} \ell_1(\lambda) &=\log L_1(\lambda;x)\\ &=\log(e^{-\lambda})+\log(\lambda^x)-\log(x!)\\ &=-\lambda+x\log\lambda-\log(x!). \end{aligned}$$

$x$ は観測済みの値なので、$\lambda$ で微分すると $-\log(x!)$ は消えます。したがってスコアは

$$\begin{aligned} U_1(\lambda) &=\frac{\partial}{\partial\lambda}\ell_1(\lambda)\\ &=-1+\frac{x}{\lambda}\\ &=\frac{x-\lambda}{\lambda}. \end{aligned}$$

最後の形では、分子が「観測した件数 $x$ と平均件数 $\lambda$ のずれ」です。観測値が想定より多ければスコアは正になり、$\lambda$ を大きくする方向を示します。

方法A：スコアの二乗の期待値で求める

定義 $I_1(\lambda)=E[U_1(\lambda)^2]$ に、上のスコアを代入します。

$$\begin{aligned} I_1(\lambda) &=E\left[\left(\frac{X-\lambda}{\lambda}\right)^2\right]\\ &=\frac{E[(X-\lambda)^2]}{\lambda^2}\\ &=\frac{\operatorname{Var}(X)}{\lambda^2}\\ &=\frac{\lambda}{\lambda^2}\\ &=\frac{1}{\lambda}. \end{aligned}$$

3行目で $E[(X-\lambda)^2]=\operatorname{Var}(X)$、4行目で Poisson分布の性質 $\operatorname{Var}(X)=\lambda$ を使いました。

方法B：対数尤度の曲率で確かめる

もう一度微分すると、

$$\frac{\partial^2}{\partial\lambda^2}\ell_1(\lambda) =-\frac{X}{\lambda^2}.$$

よって、Fisher情報量のもう一つの表現から

$$\begin{aligned} I_1(\lambda) &=-E\left[\frac{\partial^2}{\partial\lambda^2}\ell_1(\lambda)\right]\\ &=-E\left[-\frac{X}{\lambda^2}\right]\\ &=\frac{E[X]}{\lambda^2}\\ &=\frac{\lambda}{\lambda^2} =\frac1\lambda. \end{aligned}$$

2通りの計算が一致しました。後者は「尤度の山が母数方向にどれだけ鋭いか」を測っていると読めます。

$n$ 個の独立な観測ではどうなるか

$X_1,\ldots,X_n\overset{\mathrm{iid}}{\sim}\operatorname{Poisson}(\lambda)$ なら、対数尤度とスコアは

$$\begin{aligned} \ell_n(\lambda) &=\left(\sum_{i=1}^n X_i\right)\log\lambda -n\lambda-\sum_{i=1}^n\log(X_i!),\\ U_n(\lambda) &=\frac{\sum_{i=1}^n X_i-n\lambda}{\lambda}. \end{aligned}$$

したがって情報量は足し算され、

$$I_n(\lambda)=nI_1(\lambda)=\frac{n}{\lambda}.$$

例えば平均コロニー数が $\lambda=4$、独立な培養皿が $n=25$ 枚なら $I_{25}(4)=25/4=6.25$ です。Cramér—Rao下限は $1/I_{25}(4)=0.16$、標準誤差の下限は $\sqrt{0.16}=0.4$ です。実際、最尤推定量 $\bar X$ は $\operatorname{Var}(\bar X)=\lambda/n=4/25=0.16$ となり、この下限に一致します。

正規分布の平均（分散既知）

次に $X\sim N(\mu,\sigma^2)$ を考えます。ここでは $\sigma^2$ は既知で、平均 $\mu$ だけが未知 とします。例えば、測定系のばらつき $\sigma^2$ が検証済みで、ある化合物処置後の平均応答だけを推定する場面です。

確率密度関数は

$$f(x\mid\mu)=\frac{1}{\sqrt{2\pi\sigma^2}} \exp\left\{-\frac{(x-\mu)^2}{2\sigma^2}\right\}.$$

1観測から、対数尤度とスコアをつくる

対数を取ると、

$$\ell_1(\mu) =-\frac12\log(2\pi\sigma^2) -\frac{(x-\mu)^2}{2\sigma^2}.$$

第2項だけが $\mu$ を含みます。平方を微分する際には

$$\frac{\partial}{\partial\mu}(x-\mu)^2 =2(x-\mu)(-1) =-2(x-\mu)$$

であることに注意すると、

$$\begin{aligned} U_1(\mu) &=\frac{\partial}{\partial\mu}\ell_1(\mu)\\ &=-\frac{1}{2\sigma^2}\{-2(x-\mu)\}\\ &=\frac{x-\mu}{\sigma^2}. \end{aligned}$$

これは「観測値と現在の平均候補の差」を、測定誤差の大きさ $\sigma^2$ で割ったものです。同じずれでも、測定が精密で $\sigma^2$ が小さいほど、平均を動かす根拠は強くなります。

方法A：スコアの二乗の期待値で求める

$$\begin{aligned} I_1(\mu) &=E\left[\left(\frac{X-\mu}{\sigma^2}\right)^2\right]\\ &=\frac{E[(X-\mu)^2]}{\sigma^4}\\ &=\frac{\operatorname{Var}(X)}{\sigma^4}\\ &=\frac{\sigma^2}{\sigma^4}\\ &=\frac1{\sigma^2}. \end{aligned}$$

ここでも $E[(X-\mu)^2]=\operatorname{Var}(X)$ を使っています。分散が小さい測定ほど $I_1(\mu)$ は大きくなります。

方法B：対数尤度の曲率で確かめる

スコアをもう一度微分すると、観測値 $X$ を含まない定数になります。

$$\frac{\partial^2}{\partial\mu^2}\ell_1(\mu) =-\frac1{\sigma^2}.$$

したがって、

$$I_1(\mu) =-E\left[-\frac1{\sigma^2}\right] =\frac1{\sigma^2}.$$

$n$ 個の独立な測定ではどうなるか

$$\begin{aligned} \ell_n(\mu) &=-\frac{n}{2}\log(2\pi\sigma^2) -\frac1{2\sigma^2}\sum_{i=1}^n(X_i-\mu)^2,\\ U_n(\mu) &=\frac{1}{\sigma^2}\sum_{i=1}^n(X_i-\mu),\\ I_n(\mu)&=\frac{n}{\sigma^2}. \end{aligned}$$

例えば既知の標準偏差が $\sigma=2$、独立な測定が $n=16$ 回なら、$I_{16}(\mu)=16/2^2=4$ です。推定量 $\bar X$ の分散は $\sigma^2/n=4/16=0.25$、標準誤差は $0.5$ であり、情報量の逆数 $1/4=0.25$ とぴったり一致します。

発展：平均と分散がともに未知なら、情報量はどうなるか

分散そのものを $v=\sigma^2$ と置くと、1観測の対数尤度は

$$\ell_1(\mu,v) =-\frac12\log(2\pi v)-\frac{(x-\mu)^2}{2v}$$

です。2つの未知母数に対するスコアは

$$U_\mu=\frac{x-\mu}{v}, \qquad U_v=-\frac1{2v}+\frac{(x-\mu)^2}{2v^2}.$$

したがって1観測の情報行列は

$$I_1(\mu,v)= \begin{pmatrix} 1/v & 0\\ 0 & 1/(2v^2) \end{pmatrix} = \begin{pmatrix} 1/\sigma^2 & 0\\ 0 & 1/(2\sigma^4) \end{pmatrix}.$$

左上は先ほどの「平均に関する情報量」です。右下は分散に関する情報量、非対角成分の0は、この正規分布のパラメータ化では平均と分散のスコアが無相関であることを意味します。0だからといって、一般の分布で常に2つの推定が独立になるわけではありません。

モデル	1観測の情報量	$n$ 観測の情報量	$\bar X$ の分散
$\operatorname{Poisson}(\lambda)$	$1/\lambda$	$n/\lambda$	$\lambda/n$
$N(\mu,\sigma^2)$（$\sigma^2$既知）	$1/\sigma^2$	$n/\sigma^2$	$\sigma^2/n$

9. Cramér–Rao不等式

9.1 不偏推定量の分散の下限

$T$ が $\theta$ の不偏推定量なら、

$$\operatorname{Var}_\theta(T) \ge \frac{1}{I_n(\theta)}.$$

これは、正則条件下で不偏推定量が超えられない精度の限界です。

9.2 導出

ここでは、独立標本 $X=(X_1,\ldots,X_n)$ の同時密度（離散なら同時確率関数）を $f_n(x\mid\theta)$、推定量を $T=T(X)$ と書きます。示したいことは、次の不等式です。

$$\operatorname{Var}_\theta(T)\ge\frac{1}{I_n(\theta)}.$$

式だけを追う前に、意味を一言でいうとこうです。

不偏推定量は、真の母数を平均として正しく狙うため、尤度の傾き（スコア）と一定以上結び付いていなければなりません。その結び付きの強さにはCauchy—Schwarz不等式による上限があるため、推定量のばらつきは無限には小さくできません。

以下の4段階を順に確認します。

段階	行うこと	得られる式
1	不偏性を母数で微分する	$\partial E[T]/\partial\theta=1$
2	密度の微分をスコアで表す	$1=E[TU_n]$
3	スコアの平均が0と使う	$\operatorname{Cov}(T,U_n)=1$
4	Cauchy—Schwarz不等式を使う	$1\le\operatorname{Var}(T)I_n(\theta)$

段階1：不偏性を微分する

$T$ が $\theta$ の不偏推定量であることは、

$$E_\theta[T]=\theta$$

を意味します。例えば標本比率 $\bar X$ は、Bernoulli分布の反応確率 $p$ に対して $E[\bar X]=p$ なので不偏です。

両辺を $\theta$ で微分すると、右辺は $\theta$ の1階関数なので

$$\frac{\partial}{\partial\theta}E_\theta[T]=1.$$

左辺は、期待値を同時密度で書き直してから微分します。

$$\begin{aligned} \frac{\partial}{\partial\theta}E_\theta[T] &=\frac{\partial}{\partial\theta} \int T(x)f_n(x\mid\theta)\,dx\\ &=\int T(x)\frac{\partial}{\partial\theta} f_n(x\mid\theta)\,dx. \end{aligned}$$

ここで $T(x)$ は観測データから作った量であり、$\theta$ の関数ではないため、微分されるのは $f_n(x\mid\theta)$ だけです。また、上の2行目では微分と積分を交換できるという正則条件を使っています。

段階2：密度の微分を「密度 × スコア」に変える

スコアは

$$U_n(\theta) =\frac{\partial}{\partial\theta}\log f_n(X\mid\theta)$$

でした。対数の微分 $\partial\log f/\partial\theta=(\partial f/\partial\theta)/f$ を逆向きに使うと、

$$\frac{\partial}{\partial\theta}f_n(x\mid\theta) =f_n(x\mid\theta) \frac{\partial}{\partial\theta}\log f_n(x\mid\theta) =f_n(x\mid\theta)U_n(\theta)$$

です。これを段階1の積分へ代入します。

$$\begin{aligned} 1 &=\int T(x)f_n(x\mid\theta)U_n(\theta)\,dx\\ &=E_\theta[TU_n(\theta)]. \end{aligned}$$

ここまでは「不偏である」という条件だけから、推定量 $T$ とスコアの積の期待値が1になることを導きました。

段階3：なぜ共分散が1になるのか

スコアの期待値は0です。実際、

$$\begin{aligned} E_\theta[U_n(\theta)] &=\int U_n(\theta)f_n(x\mid\theta)\,dx\\ &=\int \frac{\partial}{\partial\theta}f_n(x\mid\theta)\,dx\\ &=\frac{\partial}{\partial\theta} \int f_n(x\mid\theta)\,dx\\ &=\frac{\partial}{\partial\theta}1=0. \end{aligned}$$

共分散の定義

$$\operatorname{Cov}(T,U_n) =E[TU_n]-E[T]E[U_n]$$

へ、$E[TU_n]=1$ と $E[U_n]=0$ を入れると、

$$\operatorname{Cov}(T,U_n) =1-E[T]\cdot0 =1.$$

これは、推定量が正しく母数を追うなら、尤度を増やす方向を示すスコアと必ず一定の関係を持つ、という意味です。

段階4：Cauchy—Schwarz不等式を使う

Cauchy—Schwarz不等式は、2つの確率変数の共分散について

$$\operatorname{Cov}(T,U_n)^2 \le \operatorname{Var}(T)\operatorname{Var}(U_n)$$

といいます。左辺は段階3から $1^2$ です。右辺の第2因子は、スコアの平均が0であるため

$$\operatorname{Var}(U_n) =E[U_n^2]-\{E[U_n]\}^2 =E[U_n^2] =I_n(\theta)$$

です。よって、

$$\begin{aligned} 1 &=\operatorname{Cov}(T,U_n)^2\\ &\le \operatorname{Var}(T)\operatorname{Var}(U_n)\\ &=\operatorname{Var}(T)I_n(\theta). \end{aligned}$$

両辺を正の $I_n(\theta)$ で割れば、

$$\operatorname{Var}(T)\ge\frac1{I_n(\theta)}.$$

これがCramér—Rao不等式です。情報量が大きいほど、許される最小分散は小さくなります。ただし、これは「どんな推定量でも達成できる分散」ではなく、正則条件を満たす不偏推定量が下回れない下限です。

通し例：Bernoulli分布で共分散を実際に計算する

$X_i\overset{\mathrm{iid}}{\sim}\operatorname{Bernoulli}(p)$、$S=\sum_{i=1}^nX_i$ とします。$S\sim\operatorname{Binomial}(n,p)$ なので、標本比率 $T=\hat p=S/n$ は

$$E[\hat p]=E[S]/n=np/n=p$$

を満たす不偏推定量です。対数尤度とスコアは

$$\begin{aligned} \ell_n(p)&=S\log p+(n-S)\log(1-p),\\ U_n(p)&=\frac{S}{p}-\frac{n-S}{1-p} =\frac{S-np}{p(1-p)}. \end{aligned}$$

ここで $S-np$ は $S$ から平均 $E[S]=np$ を引いたものです。したがって共分散は、定数を引いても変わらないことを使って

$$\begin{aligned} \operatorname{Cov}(\hat p,U_n) &=\operatorname{Cov}\left( \frac{S}{n},\frac{S-np}{p(1-p)} \right)\\ &=\frac{1}{np(1-p)}\operatorname{Cov}(S,S-np)\\ &=\frac{1}{np(1-p)}\operatorname{Var}(S)\\ &=\frac{np(1-p)}{np(1-p)}=1. \end{aligned}$$

一般の導出で出てきた $\operatorname{Cov}(T,U_n)=1$ が、具体例でも確かめられました。

続いて、標本比率の分散は

$$\operatorname{Var}(\hat p) =\operatorname{Var}\left(\frac1n\sum_iX_i\right) =\frac1{n^2}\sum_i\operatorname{Var}(X_i) =\frac{p(1-p)}n.$$

一方、

$$\frac1{I_n(p)} =\frac1{n/\{p(1-p)\}} =\frac{p(1-p)}n.$$

両者が一致するため、標本比率はこの正則モデルでCramér—Rao下限へ達します。

9.3 $g(\theta)$ の不偏推定

今までは $E[T]=\theta$、すなわち母数そのものを推定する場合でした。$g(\theta)$ の不偏推定量では、出発点だけが

$$E_\theta[T]=g(\theta)$$

に変わります。これを微分すると、段階1〜3と同じ計算により

$$g'(\theta) =\frac{\partial}{\partial\theta}E_\theta[T] =E_\theta[TU_n(\theta)] =\operatorname{Cov}_\theta(T,U_n(\theta))$$

です。最後の等号では、やはり $E[U_n]=0$ を使いました。Cauchy—Schwarz不等式を適用すると、

$$\{g'(\theta)\}^2 \le\operatorname{Var}(T)I_n(\theta),$$

よって、

$$\operatorname{Var}(T) \ge\frac{\{g'(\theta)\}^2}{I_n(\theta)}.$$

$g'(\theta)$ は「推定したい量が、母数の変化へどれだけ敏感か」です。例えば $g(\theta)=\log\theta$ なら $g'(\theta)=1/\theta$ であり、母数の尺度を変えると分散の下限も変わります。単に $1/I_n$ を機械的に使わず、何を推定したいかが $g(\theta)$ かを確認します。

9.4 等号成立と有効推定量

Cramér—Rao不等式で使ったCauchy—Schwarz不等式は、実際には中心化した2変数

$$A=T-E[T]=T-\theta, \qquad B=U_n(\theta)-E[U_n(\theta)]=U_n(\theta)$$

へ適用していました。不偏性 $E[T]=\theta$ と、スコアの性質 $E[U_n]=0$ を使ったため、この形になります。

では、なぜ等号のときに $A$ と $B$ が比例するのでしょうか。任意の定数 $c$ について、平方の期待値は必ず非負なので、

$$\begin{aligned} 0 &\le E[(A-cB)^2]\\ &=E[A^2]-2cE[AB]+c^2E[B^2]\\ &=\operatorname{Var}(A)-2c\operatorname{Cov}(A,B) +c^2\operatorname{Var}(B). \end{aligned}$$

ここで、右辺をもっとも小さくする

$$c=\frac{\operatorname{Cov}(A,B)}{\operatorname{Var}(B)}$$

を選ぶと、

$$0\le \operatorname{Var}(A) -\frac{\operatorname{Cov}(A,B)^2}{\operatorname{Var}(B)}.$$

両辺へ $\operatorname{Var}(B)$ を掛ければCauchy—Schwarz不等式です。そして等号が成り立つのは、最初の平方の期待値が0、すなわち

$$E[(A-cB)^2]=0$$

となるときだけです。非負な確率変数の期待値が0なら、その確率変数は確率1で0です。よって

$$A=cB \quad\text{（確率1で）}$$

が等号成立の必要十分条件です。

これをCramér—Rao不等式へ戻します。ここでは

$$\operatorname{Cov}(A,B) =\operatorname{Cov}(T,U_n)=1, \qquad \operatorname{Var}(B) =\operatorname{Var}(U_n)=I_n(\theta)$$

なので、比例定数は

$$c=\frac{1}{I_n(\theta)}$$

と決まります。したがって、下限へ達する条件は単に「ある定数倍」ではなく、

$$T-\theta =\frac{U_n(\theta)}{I_n(\theta)} \quad\text{（確率1で）}$$

と書けることです。

これは「推定誤差 $T-\theta$ がスコアとまったく同じ方向にだけ揺れる」ことを意味します。一般には成り立たない強い条件です。

Poisson分布では $\bar X$ は $\lambda$ の不偏推定量で、

$$\operatorname{Var}(\bar X)=\frac{\lambda}{n}.$$

一方、

$$\frac1{I_n(\lambda)} =\frac1{n/\lambda} =\frac{\lambda}{n}.$$

実際、$U_n(\lambda)=n(\bar X-\lambda)/\lambda$ なので、

$$\bar X-\lambda=\frac{\lambda}{n}U_n(\lambda)$$

一方、$I_n(\lambda)=n/\lambda$ なので

$$\frac{U_n(\lambda)}{I_n(\lambda)} =\frac{\lambda}{n}U_n(\lambda) =\bar X-\lambda.$$

これは上の等号条件の形そのものです。したがって分散も下限と一致し、$\bar X$ は有効推定量です。

$g(\theta)$ を推定する場合の等号条件

$E[T]=g(\theta)$ のときは $\operatorname{Cov}(T,U_n)=g'(\theta)$ でした。したがって下限へ達する条件は、$T-g(\theta)=\{g'(\theta)/I_n(\theta)\}U_n(\theta)$ です。母数そのものを推定する場合は $g'(\theta)=1$ となり、本文の式に戻ります。

Rao–Blackwell化とLehmann–Scheffé定理

不偏推定量 $U$ と十分統計量 $T$ があるとき、

$$\tilde U=E[U\mid T]$$

も不偏で、全分散公式から

$$\operatorname{Var}(\tilde U)\le\operatorname{Var}(U)$$

です。さらに $T$ が完全十分統計量なら、$T$ の関数として書ける不偏推定量は一意な一様最小分散不偏推定量です。

10. 一致性、漸近正規性、漸近有効性

10.1 最尤推定量の一致性

最尤推定量の性質は、次の順番で捉えると混乱しにくくなります。

性質	日常語での意味	この節で使う主な道具
一致性	標本数を増やすと真値へ近づく	大数の法則
漸近正規性	真値の近くで、誤差が正規分布らしく揺れる	中心極限定理、Taylor展開
漸近有効性	大標本では理論上もっとも小さい分散に近づく	Fisher情報量

まず一致性です。対数尤度を標本数 $n$ で割った

$$\frac1n\ell_n(\theta)$$

は「1観測当たり、候補 $\theta$ がどれだけデータを説明できるか」を表します。$n$ が大きくなると、偶然による上下が平均化されます。

大数の法則により、1標本当たりの対数尤度は

$$\frac1n\ell_n(\theta) =\frac1n\sum_{i=1}^n\log f(X_i\mid\theta) \xrightarrow{p} E_{\theta_0}[\log f(X\mid\theta)].$$

右辺が真値 $\theta_0$ で一意に最大なら、対数尤度の最大点も真値へ近づき、

$$\hat\theta_{\mathrm{ML}}\xrightarrow{p}\theta_0$$

となります。厳密には一様収束、識別可能性、母数空間などの条件が必要です。

ここで識別可能性とは、異なる母数が同じ分布を作らないことです。たとえばPoisson分布なら平均が異なれば分布も異なるため、$\lambda$ は識別できます。初読では「大標本では偶然の揺れより平均的な当てはまりが勝つ」と理解すれば十分です。

10.2 漸近正規性をTaylor展開から導く

この節で最終的に知りたいのは、推定値の誤差がおおよそどのくらいかです。結論だけ先に書くと、

$$\hat\theta_{\mathrm{ML}} \mathrel{\dot\sim} N\left(\theta_0,\frac1{nI_1(\theta_0)}\right).$$

右辺の分散は $1/n$ の速さで小さくなります。導出では「最尤推定量ではスコアが0」という事実を、真値 $\theta_0$ の近くでTaylor展開します。

最尤推定量はスコア方程式

$$U_n(\hat\theta)=0$$

を満たすとします。真値 $\theta_0$ の周りで1次Taylor展開すると、

$$0 =U_n(\theta_0) +(\hat\theta-\theta_0)U_n'(\tilde\theta)$$

となる $\tilde\theta$ が $\theta_0$ と $\hat\theta$ の間に存在します。整理して、

$$\sqrt n(\hat\theta-\theta_0) = \frac{U_n(\theta_0)/\sqrt n} {-U_n'(\tilde\theta)/n}.$$

この式は「推定誤差 = 真値で残った傾き ÷ 尤度の曲がり具合」と読めます。分子と分母をそれぞれ別の定理で扱います。

分子：真値でのスコアの揺れ

スコアは独立な1標本スコアの和です。

$$U_n(\theta_0)=\sum_{i=1}^nU_1^{(i)}(\theta_0).$$

各項は平均0、分散 $I_1(\theta_0)$ を持つので、中心極限定理から

$$\frac{U_n(\theta_0)}{\sqrt n} \xrightarrow{d}N\{0,I_1(\theta_0)\}.$$

ここで $\sqrt n$ で割るのは、独立な和の標準偏差が $\sqrt n$ の大きさになるためです。

分母：尤度の曲がり具合

$U_n'(\theta)$ も1標本ごとの2階微分の和です。$\hat\theta$ が一致して $\theta_0$ へ近づくなら、その間にある $\tilde\theta$ も $\theta_0$ へ近づきます。大数の法則により、

$$-\frac1nU_n'(\tilde\theta) \xrightarrow{p}I_1(\theta_0).$$

つまり、分母はランダムに見えても大標本では Fisher情報量という正の定数へ落ち着きます。

最後に分子と分母を合わせる

分子は正規分布へ、分母は定数へ近づくため、Slutskyの定理から

$$\sqrt n(\hat\theta-\theta_0) \xrightarrow{d} N\left(0,\frac1{I_1(\theta_0)}\right).$$

したがって大標本では、

$$\hat\theta \mathrel{\dot\sim} N\left( \theta_0,\frac1{nI_1(\theta_0)} \right).$$

Poisson分布で抽象式を確かめる

Poisson分布では $\hat\lambda_{\mathrm{ML}}=\bar X$、$I_1(\lambda)=1/\lambda$ です。したがって一般式へ代入すると、

$$\sqrt n(\bar X-\lambda) \xrightarrow{d}N(0,\lambda).$$

これは中心極限定理

$$\sqrt n(\bar X-E[X]) \xrightarrow{d}N\{0,\operatorname{Var}(X)\}$$

へ $E[X]=\operatorname{Var}(X)=\lambda$ を代入した結果と同じです。最尤推定の一般理論が、既知の標本平均の理論へ戻ることを確認できます。

10.3 漸近分散と漸近有効性

$$\sqrt n(\hat\theta-\theta) \xrightarrow{d}N\{0,V(\theta)\}$$

の $V(\theta)$ を漸近分散と呼びます。正則な推定問題では、漸近分散が $I_1(\theta)^{-1}$ に達する推定量を漸近有効といいます。

ここで混同しやすい点は、$V(\theta)$ 自体は $\sqrt n$ を掛けた誤差の分散だということです。元の推定量の大標本での分散は、

$$\operatorname{Var}(\hat\theta) \approx\frac{V(\theta)}{n}$$

と読みます。最尤推定量なら $V(\theta)=1/I_1(\theta)$ なので、標準誤差はおおよそ $1/\sqrt{nI_1(\theta)}$ です。

モデル	$I_1(\theta)$	最尤推定量の近似分散
$\operatorname{Poisson}(\lambda)$	$1/\lambda$	$\lambda/n$
$N(\mu,\sigma^2)$（$\sigma^2$既知）	$1/\sigma^2$	$\sigma^2/n$

最尤推定量は適切な正則条件のもとで、一致性、漸近正規性、漸近有効性を持ちます。ただし小標本、境界母数、混合モデル、識別不能、強い外れ値では近似が悪いことがあります。

M推定量とサンドイッチ分散

推定方程式

$$\sum_{i=1}^n\psi(X_i,\theta)=0$$

で定まるM推定量は、最尤推定量を含む広いクラスです。モデルが完全には正しくなくても、

$$\sqrt n(\hat\theta-\theta_0) \xrightarrow{d} N(0,A^{-1}BA^{-\mathsf T})$$

となることがあります。ここで

$$A=E\left[-\frac{\partial\psi}{\partial\theta^\mathsf T}\right], \qquad B=E[\psi\psi^\mathsf T].$$

1母数なら何に当たるか

行列の逆は数の逆数になり、サンドイッチ分散は $B/A^2$ です。最尤法でモデルが正しければ $A=B=I_1(\theta)$ となり、$B/A^2=1/I_1(\theta)$ に戻ります。

11. 発展的なBayes推定

経験Bayes

事前分布の超母数を外部から固定せず、多数の群のデータから推定してから各群を更新する方法です。

例えば複数遺伝子の発現差や複数施設の有害事象率では、群ごとの生の推定値 $\hat\theta_j$ を、全体平均へ適度に縮小できます。小標本群ほど強く縮み、大標本群はデータを保ちます。

利点は計算の軽さです。一方、超母数を推定した不確実性を無視しやすく、群数が少ないと不安定です。

階層Bayes

階層モデルでは、超母数も未知量として事前分布を置きます。

$$\begin{aligned} Y_{ij}\mid\theta_j,\sigma^2&\sim N(\theta_j,\sigma^2),\\ \theta_j\mid\mu,\tau^2&\sim N(\mu,\tau^2),\\ (\mu,\tau)&\sim \pi(\mu,\tau). \end{aligned}$$

観測値、群別母数、集団母数という3層を同時に推定します。個体差や施設差を持つ薬物動態、動物ごとの細胞データ、プレート差を含むassayに向きます。

部分プーリングを数値で見る

群 $j$ の観測平均を $\bar Y_j$、その分散を $v_j$ とし、

$$\bar Y_j\mid\theta_j\sim N(\theta_j,v_j), \qquad \theta_j\sim N(\mu,\tau^2)$$

とします。正規分布どうしを掛けて平方完成すると、事後平均は

$$E[\theta_j\mid\bar Y_j] =\frac{\tau^2}{\tau^2+v_j}\bar Y_j +\frac{v_j}{\tau^2+v_j}\mu$$

になります。例えば全体平均 $\mu=10$、群間分散 $\tau^2=4$、ある群の観測平均 $\bar Y_j=16$、観測分散 $v_j=9$ なら、

$$E[\theta_j\mid\bar Y_j] =\frac4{13}\times16+\frac9{13}\times10 =\frac{154}{13} \approx11.85.$$

生の平均16をそのまま使わず、情報の少ない群を全体平均10へ縮めています。$v_j$ が小さく、群内データが精密になるほど $\bar Y_j$ の重みが大きくなります。

正規分布の積から重みを導く

事後分布の指数部は $-(\bar Y_j-\theta_j)^2/(2v_j)-(\theta_j-\mu)^2/(2\tau^2)$ です。$\theta_j$ の2次式としてまとめると、事後精度は $1/v_j+1/\tau^2$、事後平均は $\{(1/v_j)\bar Y_j+(1/\tau^2)\mu\}/(1/v_j+1/\tau^2)$ となります。分母分子へ $v_j\tau^2$ を掛けると本文の加重平均になります。

客観BayesとJeffreys事前分布

主観的な事前情報を弱める方法の1つが、

$$\pi_J(\theta)\propto\sqrt{I_1(\theta)}$$

で定義されるJeffreys事前分布です。滑らかな1対1変換に対して不変です。

ただし積分が1にならない不適切事前分布になることがあります。事後分布が適切か、境界で発散しないかを確認しなければなりません。「客観」は仮定がないという意味ではありません。

Bernstein–von Misesの見方

正則条件と大標本のもとで、事後分布は最尤推定量を中心とする正規分布へ近づき、

$$\pi(\theta\mid X) \approx N\left(\hat\theta_{\mathrm{ML}},\{nI_1(\theta_0)\}^{-1}\right).$$

頻度論とBayesの区間が近づく理由ですが、高次元、境界、混合分布、非識別モデルでは成立しないことがあります。

12. 多次元への拡張

12.1 スコアベクトルとFisher情報行列

$k$ 次元母数 $\theta=(\theta_1,\ldots,\theta_k)^\mathsf T$ では、

$$U(\theta) =\nabla_\theta\ell(\theta) = \begin{pmatrix} \partial\ell/\partial\theta_1\\ \vdots\\ \partial\ell/\partial\theta_k \end{pmatrix}.$$

Fisher情報行列は、

$$I(\theta) =E[U(\theta)U(\theta)^\mathsf T] =-E[\nabla_\theta^2\ell(\theta)].$$

対角成分は各母数自身の情報量、非対角成分は母数間の推定上の結び付きを表します。

たとえば、EC50とHill係数を同時に推定するとき、情報行列の左上・右下はそれぞれの推定精度に関係し、非対角成分は「EC50を変えた影響をHill係数でどの程度補えてしまうか」を表します。非対角成分が大きいほど、2つを別々に決めにくくなります。

2母数なら各要素をどう計算するか

$U=(U_1,U_2)^\mathsf T$ なら、$I_{11}=E[U_1^2]$、$I_{22}=E[U_2^2]$、$I_{12}=I_{21}=E[U_1U_2]$ です。1母数では1行1列になり、普通の情報量 $E[U^2]$ に戻ります。

2×2行列を実際に逆行列へする

例えば、2母数の情報行列が

$$I= \begin{pmatrix} 4&1\\ 1&2 \end{pmatrix}$$

だったとします。2×2行列

$$\begin{pmatrix}a&b\\c&d\end{pmatrix}^{-1} =\frac1{ad-bc} \begin{pmatrix}d&-b\\-c&a\end{pmatrix}$$

を使うと、行列式は $4\times2-1\times1=7$ なので、

$$I^{-1} =\frac17 \begin{pmatrix} 2&-1\\ -1&4 \end{pmatrix}.$$

1母数なら「情報量4の逆数は分散下限 $1/4$」でした。2母数では逆行列の対角成分 $2/7$ と $4/7$ が各推定量の分散下限に対応し、非対角成分 $-1/7$ が2つの推定誤差の結び付きを表します。

12.2 多次元Cramér–Rao不等式

不偏推定量ベクトル $\hat\theta$ について、

$$\operatorname{Cov}(\hat\theta)-I_n(\theta)^{-1}\succeq0$$

です。$\succeq0$ は左辺が半正定値、すなわち任意のベクトル $a$ に対して

$$a^\mathsf T\operatorname{Cov}(\hat\theta)a \ge a^\mathsf TI_n(\theta)^{-1}a$$

であることを意味します。

左辺の $a^\mathsf T\operatorname{Cov}(\hat\theta)a$ は、線形結合 $a^\mathsf T\hat\theta$ の分散です。つまり多次元版は「どの母数の組合せを見ても、その推定誤差の分散は情報行列の逆数で決まる下限を下回れない」と言っています。1母数の不等式を、あらゆる方向へ同時に拡張したものです。

単変量なら何に当たるか

$a$ も行列も数になるため、$\operatorname{Var}(\hat\theta)\ge1/I_n(\theta)$ です。ベクトル $a$ は「複数母数のどの線形結合を見たいか」を選ぶ役割です。

12.3 多次元最尤推定量の漸近正規性

$$\sqrt n(\hat\theta_{\mathrm{ML}}-\theta_0) \xrightarrow{d} N_k\{0,I_1(\theta_0)^{-1}\}.$$

ここで共分散行列の対角成分は各推定量の漸近分散、非対角成分は推定量どうしの漸近共分散です。

この式の役割は、1母数の $1/\{nI_1(\theta)\}$ を行列版の $I_1(\theta_0)^{-1}/n$ に置き換えることです。実務ではこの逆行列の対角成分から各母数の標準誤差を、非対角成分から推定値間のトレードオフを読みます。

式を単変量へ戻して確認する

$k=1$ なら、$N_k$ は通常の正規分布 $N$、共分散行列 $I_1^{-1}$ は数 $1/I_1$ です。したがって $\sqrt n(\hat\theta-\theta_0)\to N(0,1/I_1)$ となり、10.2節の式へそのまま戻ります。

12.4 邪魔母数とプロファイル尤度

関心母数を $\psi$、それ以外を $\lambda$ と分けます。各 $\psi$ に対して $\lambda$ を最大化した

$$L_p(\psi)=L\{\psi,\hat\lambda(\psi)\}$$

をプロファイル尤度といいます。非線形薬効モデルでEC50に関心があり、最大反応やHill係数が邪魔母数になる場合に使えます。

操作は3段階です。1. EC50候補 $\psi$ を1つ固定する、2. その条件で最大反応・Hill係数などを最も合うように調整する、3. 残った尤度をEC50候補どうしで比べる、です。ほかの母数を無視する方法ではなく、各候補で最も有利な条件を与えたうえでEC50を比べる方法です。

Schur補完で見る邪魔母数の情報損失

情報行列を

$$I(\theta)= \begin{pmatrix} I_{\psi\psi}&I_{\psi\lambda}\\ I_{\lambda\psi}&I_{\lambda\lambda} \end{pmatrix}$$

と分けると、邪魔母数未知のときに $\psi$ へ残る有効情報は

$$I_{\psi\psi\cdot\lambda} =I_{\psi\psi} -I_{\psi\lambda}I_{\lambda\lambda}^{-1}I_{\lambda\psi}.$$

第2項だけ情報が失われます。母数間の相関が強いほど、関心母数の推定が不安定になります。

2母数の数の計算に戻す

各ブロックが数なら $I_{\psi\psi\cdot\lambda}=I_{\psi\psi}-I_{\psi\lambda}^2/I_{\lambda\lambda}$ です。非対角成分が0なら情報損失はありません。

13. 医薬・生命科学・情報科学での読み方

事例	データ	推定対象	注意点
GPCR濃度反応	濃度ごとの応答	EC50、最大反応、Hill係数	非線形性、境界、プレート差
薬物動態	時点別濃度	CL、V、個体間分散	反復測定、階層構造、BLQ
有害事象	曝露時間と件数	発生率	過分散、追跡時間の差
バイオマーカー陽性率	陽性数と総数	真の陽性率	小標本、施設差、検査誤差
RNA-seq	遺伝子別カウント	平均発現、分散、効果量	負の二項、ライブラリサイズ、多重性
scRNA-seq	細胞×遺伝子	細胞状態の効果	個体が実験単位、擬似反復
機械学習	訓練データと損失	重み、予測確率	汎化、分布外、校正

13.1 MAP推定と正則化

事後最大化は

$$\hat\theta_{\mathrm{MAP}} =\arg\max_\theta\{\log L(\theta)+\log\pi(\theta)\}.$$

負号を付ければ、負の対数尤度と罰則項の最小化です。正規事前分布はL2正則化、Laplace事前分布はL1正則化に対応します。

ただしニューラルネットワークの重みについてMAPを求めただけでは、Bayes的な予測不確実性を積分したことにはなりません。

13.2 EC50推定で確認すること

1. 独立な実験単位は何か
2. 応答の分散は濃度で一定か
3. 最大反応や下限を固定する根拠があるか
4. EC50とHill係数が強く相関していないか
5. 点推定だけでなくプロファイル尤度やbootstrap区間を示したか

未発表研究データを公開記事へ載せる場合は、具体的な化合物名、数値、図、共同研究情報を一般化してから使用してください。

14. 推定法の選び方

状況	第一候補	理由
手計算で初期値が必要	モーメント法	簡単で高速
標準的な正則モデル・十分な標本	最尤法	一致性と漸近効率が期待できる
小標本で外部情報がある	Bayes法	事前情報と母数不確実性を統合
多数の似た群を同時推定	経験Bayes・階層Bayes	部分プーリングでMSEを下げる
不均一な精度の不偏推定値を統合	逆分散重みBLUE	線形不偏クラスで分散最小
外れ値やモデル誤指定が懸念	ロバストM推定・サンドイッチ分散	尤度仮定への依存を緩和

15. 答案で使える確認手順

1. 分布と母数空間、支持範囲を書く
2. 求める対象が母数、推定量、推定値のどれか確認する
3. 十分統計量なら因子分解を明示する
4. 最尤法なら対数尤度、1階条件、最大確認を書く
5. Bayes法なら事前×尤度の指数を整理し、分布族を同定する
6. 不偏性は期待値、MSEは分散とバイアスに分ける
7. Cramér–Raoでは正則条件と情報量が1標本か全標本かを書く
8. 漸近分布では $\sqrt n$ で尺度化し、使った定理を書く
9. 多次元式では行列の次元と、対角・非対角の意味を説明する

16. 数理統計の問題

問題の型	最初に書くもの	最後に確認するもの
十分統計量	同時確率・同時密度	母数を含む部分が統計量だけの関数か
最尤推定	尤度と許される母数の範囲	微分だけでなく最大になる範囲か
Bayes更新	事前分布 × 尤度	指数と母数化（尺度・率）が合うか
分散・MSE	$E[T]$ と $\operatorname{Var}(T)$	$\operatorname{MSE}=\operatorname{Var}+\operatorname{Bias}^2$
漸近理論	$\sqrt n(\hat\theta-\theta)$	CLT・LLN・Slutskyの役割

問題1：Bernoulli標本の十分統計量

$X_1,\ldots,X_n\overset{\mathrm{iid}}{\sim}\operatorname{Bernoulli}(p)$ とする。$T=\sum_iX_i$ が $p$ の十分統計量であることを因子分解定理で示せ。

解答1

方針： 因子分解定理では、同時確率関数を「母数 $p$ と統計量 $T$ の関数」と「データだけの関数」の積に分けます。

同時確率関数は

$$\prod_i p^{x_i}(1-p)^{1-x_i} =p^{\sum_i x_i}(1-p)^{n-\sum_i x_i}.$$

$p$ を含む部分は $T=\sum_i x_i$ だけの関数で、残りは $h(x)=1$ と置けます。因子分解定理より $T$ は十分です。

答案で気をつけること： 「和だけに依存する」と書くだけでなく、$g_p(T)h(x)$ の形を示します。

問題2：一様分布の十分統計量

$X_1,\ldots,X_n\overset{\mathrm{iid}}{\sim}U(0,\theta)$ とする。十分統計量を1つ求め、理由を示せ。

解答2

方針： 一様分布では、各観測値が $\theta$ 未満であるという条件を、最大値1つの条件へまとめます。

$$L(\theta) =\theta^{-n}\prod_i I(0<x_i<\theta) =\theta^{-n}I(0<x_{(n)}<\theta).$$

2行目では、すべての $x_i$ が $\theta$ 未満であることと、最大値 $x_{(n)}$ が $\theta$ 未満であることが同値であることを使いました。したがって

$$L(\theta) =\underbrace{\theta^{-n}I(x_{(n)}<\theta)}_{g_\theta\{x_{(n)}\}} \underbrace{I(0<x_{(1)})}_{h(x)}$$

と因子分解でき、$X_{(n)}$ が十分統計量です。

答案で気をつけること： 最大値を答えるだけでなく、支持範囲の積が1個の指示関数へまとまる途中を示します。

問題3：指数型分布族と自然母数

Poisson分布とBernoulli分布を指数型分布族の形に直し、自然母数と1標本の十分統計量を答えよ。

解答3

方針： 指数型分布族の標準形 $h(x)\exp\{\eta T(x)-A(\eta)\}$ と見比べ、$x$ に掛かる量を自然母数 $\eta$ と読み取ります。

Poisson分布は

$$f(x\mid\lambda) =\frac1{x!}\exp\{x\log\lambda-\lambda\} =\frac1{x!}\exp\{x\eta-e^\eta\}$$

なので $\eta=\log\lambda$、$T(x)=x$、$A(\eta)=e^\eta$ です。Bernoulli分布は

$$f(x\mid p) =\exp\left\{x\log\frac p{1-p}+\log(1-p)\right\} =\exp\{x\eta-\log(1+e^\eta)\}$$

なので $\eta=\log\{p/(1-p)\}$、$T(x)=x$、$A(\eta)=\log(1+e^\eta)$ です。2つの分布とも、標本全体では $\sum_iX_i$ が十分統計量になります。

答案で気をつけること： 自然母数は元の母数と同じとは限りません。logやlogitへの変換を明示します。

問題4：期待値母数

Bernoulli分布の自然母数を $\eta$ とする。$A(\eta)=\log(1+e^\eta)$ を微分し、期待値母数が $p$ になることを示せ。

解答4

方針： 期待値母数は $A'(\eta)$ です。まず普通に微分し、次に $e^\eta$ を $p$ で置き換えます。

$$\begin{aligned} A'(\eta) &=\frac{\partial}{\partial\eta}\log(1+e^\eta)\\ &=\frac{1}{1+e^\eta}\cdot e^\eta\\ &=\frac{e^\eta}{1+e^\eta}. \end{aligned}$$

$\eta=\log\{p/(1-p)\}$ なので $e^\eta=p/(1-p)$ です。したがって、

$$A'(\eta) =\frac{p/(1-p)}{1+p/(1-p)} =p.$$

さらに $A'(\eta)=p$ をもう一度 $\eta$ で微分すると、

$$A''(\eta) =\frac{e^\eta}{(1+e^\eta)^2} =p(1-p)$$

で、Bernoulli分布の分散に一致します。

答案で気をつけること： 微分後に自然母数から元の母数へ戻します。

問題5：Gamma分布のモーメント推定

形状 $\alpha$、尺度 $\beta$ のGamma分布から標本を得た。$\bar X$ と $m_2=n^{-1}\sum_i(X_i-\bar X)^2$ を用いてモーメント推定量を求めよ。

解答5

方針： Gamma分布の平均と分散を標本の平均・分散へ等置し、2本の式を連立して $\alpha,\beta$ を解きます。ここで $\beta$ は尺度です。

$$\bar X=\alpha\beta,\qquad m_2=\alpha\beta^2$$

と置きます。第2式を第1式で割ると、

$$\hat\beta=\frac{m_2}{\bar X}.$$

これを第1式へ戻して、

$$\hat\alpha=\frac{\bar X}{\hat\beta} =\frac{\bar X^2}{m_2}.$$

なお $m_2=n^{-1}\sum_i(X_i-\bar X)^2$ はモーメント法で使う標本2次中心モーメントです。不偏分散の分母 $n-1$ とは別物です。

答案で気をつけること： Gamma分布の第2母数が尺度か率かを先に確認します。

問題6：一様分布の最尤推定

$X_i\sim U(0,\theta)$ とする。$\theta$ の最尤推定量を求め、一致性を分布関数から示せ。

解答6

方針： 尤度が正になる母数範囲を先に決め、その範囲で $\theta^{-n}$ がどちらへ動くかを見ます。微分は不要です。

尤度は $\theta^{-n}I(\theta\ge X_{(n)})$ なので、

$$\hat\theta_{\mathrm{ML}}=X_{(n)}.$$

$\theta<X_{(n)}$ では尤度は0です。一方、$\theta\ge X_{(n)}$ では $\theta^{-n}$ は $\theta$ が大きくなるほど小さくなります。したがって尤度が正である範囲の左端 $X_{(n)}$ で最大になります。

$0<\varepsilon<\theta$ に対して、

$$\begin{aligned} P(|X_{(n)}-\theta|>\varepsilon) &=P(X_{(n)}<\theta-\varepsilon)\\ &=\left(\frac{\theta-\varepsilon}{\theta}\right)^n \to0. \end{aligned}$$

よって $X_{(n)}\xrightarrow{p}\theta$ です。

答案で気をつけること： $X_{(n)}\le\theta$ が常に成り立つため、絶対値の事象が片側だけになることを説明します。

問題7：一様分布の推定量とMSE

$X_i\sim U(0,\theta)$ とする。$X_{(n)}$ と $\tilde\theta=(n+1)X_{(n)}/n$ のバイアスとMSEを求めよ。

解答7

方針： まず最大値の平均・分散からバイアスを出し、$\operatorname{MSE}=\operatorname{Var}+\operatorname{Bias}^2$ をそのまま使います。

最大値について、

$$E[X_{(n)}]=\frac{n}{n+1}\theta,\qquad \operatorname{Var}(X_{(n)}) =\frac{n\theta^2}{(n+1)^2(n+2)}.$$

したがって、

$$\operatorname{Bias}(X_{(n)})=-\frac{\theta}{n+1}$$

で、

$$\begin{aligned} \operatorname{MSE}(X_{(n)}) &=\frac{n\theta^2}{(n+1)^2(n+2)} +\left(-\frac{\theta}{n+1}\right)^2\\ &=\frac{n\theta^2+(n+2)\theta^2} {(n+1)^2(n+2)}\\ &=\frac{2\theta^2}{(n+1)(n+2)}. \end{aligned}$$

$\tilde\theta$ は不偏で、

$$\operatorname{MSE}(\tilde\theta) =\operatorname{Var}(\tilde\theta) =\left(\frac{n+1}{n}\right)^2 \frac{n\theta^2}{(n+1)^2(n+2)} =\frac{\theta^2}{n(n+2)}.$$

答案で気をつけること： 不偏推定量のMSEだけが分散と一致します。$X_{(n)}$ ではバイアス平方を足します。

問題8：多項分布の最尤推定

$(X_1,\ldots,X_k)\sim\operatorname{Multinomial}(n;p_1,\ldots,p_k)$ とする。$\sum_jp_j=1$ のもとで最尤推定量を求めよ。

解答8

方針： $p_1,\ldots,p_k$ は和が1という制約を持つため、Lagrange未定乗数法を使います。以下の $r$ は制約に付ける未定乗数で、確率母数とは別の記号です。

定数を省いた対数尤度は $\ell=\sum_jx_j\log p_j$ です。Lagrange関数

$$Q=\sum_jx_j\log p_j+r\left(1-\sum_jp_j\right)$$

を作ると、

$$\frac{\partial Q}{\partial p_j} =\frac{x_j}{p_j}-r=0$$

より $p_j=x_j/r$ です。制約 $\sum_jp_j=1$ へ代入すると、

$$1=\sum_j\frac{x_j}{r} =\frac{\sum_jx_j}{r} =\frac nr$$

なので $r=n$、したがって

$$\hat p_j=\frac{x_j}{n}.$$

答案で気をつけること： $p_k=1-\sum_{j<k}p_j$ という制約を無視して独立に微分してはいけません。

問題9：Beta–Binomialの事後分布

$X\mid p\sim\operatorname{Binomial}(n,p)$、$p\sim\operatorname{Beta}(\alpha,\beta)$ とする。事後分布と二乗誤差損失でのBayes推定量を求めよ。

解答9

方針： 事後分布は「尤度 × 事前分布」です。$p$ を含まない二項係数やBeta関数の定数は比例記号の中へ入れます。

尤度と事前分布を $p$ のべきとして書くと、

$$p(x\mid p)\propto p^x(1-p)^{n-x}, \qquad \pi(p)\propto p^{\alpha-1}(1-p)^{\beta-1}.$$

両者を掛けると、

$$\pi(p\mid x) \propto p^{x+\alpha-1}(1-p)^{n-x+\beta-1}.$$

したがって、

$$p\mid x\sim\operatorname{Beta}(\alpha+x,\beta+n-x).$$

二乗誤差損失では事後平均が最適なので、

$$\hat p_B=\frac{\alpha+x}{\alpha+\beta+n}.$$

答案で気をつけること： MAPではなく事後平均です。損失関数で答えが変わります。

問題10：Gamma–Poissonの事後分布

$X_i\mid\lambda\sim\operatorname{Poisson}(\lambda)$、$\lambda\sim\operatorname{Gamma}(\alpha,\beta)$ とする。第2母数を率として、事後分布を求めよ。

解答10

方針： Gamma分布を「形状 $a$、率 $b$」で $f(\lambda)\propto\lambda^{a-1}e^{-b\lambda}$ と書き、尤度と掛けた後の指数を読み取ります。

尤度は $\lambda^{\sum_i x_i}e^{-n\lambda}$ に比例し、事前分布は $\lambda^{\alpha-1}e^{-\beta\lambda}$ に比例します。掛け合わせると、

$$\begin{aligned} \pi(\lambda\mid x) &\propto \lambda^{\sum_i x_i}e^{-n\lambda} \lambda^{\alpha-1}e^{-\beta\lambda}\\ &=\lambda^{\alpha+\sum_i x_i-1} \exp\{-(\beta+n)\lambda\}. \end{aligned}$$

これは形状が $\alpha+\sum_i x_i$、率が $\beta+n$ のGamma分布なので、

$$\lambda\mid x \sim\operatorname{Gamma} \left(\alpha+\sum_i x_i,\beta+n\right).$$

答案で気をつけること： 曝露時間が観測ごとに違うなら $n$ ではなく総曝露時間が率母数へ加わります。

問題11：逆分散重みBLUE

独立な不偏推定量 $Y_1,Y_2$ の分散がそれぞれ $\sigma_1^2,\sigma_2^2$ である。$aY_1+(1-a)Y_2$ の分散を最小にする $a$ を求めよ。

解答11

方針： 重みの和を1にして不偏性を保ち、重み $a$ について分散を最小化します。

不偏性は、重みの和が $a+(1-a)=1$ なので自動的に保たれます。独立性より分散は、

$$V(a)=a^2\sigma_1^2+(1-a)^2\sigma_2^2.$$

微分して、

$$V'(a)=2a\sigma_1^2-2(1-a)\sigma_2^2=0.$$

よって、

$$\begin{aligned} a\sigma_1^2&=(1-a)\sigma_2^2\\ a(\sigma_1^2+\sigma_2^2)&=\sigma_2^2\\ a&=\frac{\sigma_2^2}{\sigma_1^2+\sigma_2^2}. \end{aligned}$$

さらに $V''(a)=2(\sigma_1^2+\sigma_2^2)>0$ なので、これは最小値です。このときの最小分散は

$$V(a) =\frac{\sigma_1^2\sigma_2^2}{\sigma_1^2+\sigma_2^2} =\frac{1}{\sigma_1^{-2}+\sigma_2^{-2}}.$$

この形は「$Y_1$ の重みに相手 $Y_2$ の分散が付く」ように見えますが、分子分母に $1/(\sigma_1^2\sigma_2^2)$ を掛ければ、

$$a=\frac{\sigma_2^2}{\sigma_1^2+\sigma_2^2} =\frac{\sigma_1^{-2}}{\sigma_1^{-2}+\sigma_2^{-2}}.$$

つまり、各測定の重みは自分自身の分散の逆数に比例します。

答案で気をつけること： 最初に独立性を使って分散の共分散項を0にしたことと、重みの和が1なので不偏性が保たれることを書きます。

問題12：Poisson分布の情報量とCR下限

$X_1,\ldots,X_n\sim\operatorname{Poisson}(\lambda)$ とする。$\lambda$ のFisher情報量と不偏推定量の分散下限を求め、$\bar X$ が限界へ達するか確認せよ。

解答12

方針： 対数尤度を微分してスコアを出し、スコアの二乗の期待値から1観測の情報量を求めます。その後に $n$ 倍し、CR下限と標本平均の分散を比べます。

1観測の対数尤度は

$$\ell_1(\lambda)=X\log\lambda-\lambda-\log(X!)$$

なので、スコアは

$$U_1(\lambda) =\frac{\partial\ell_1}{\partial\lambda} =\frac{X}{\lambda}-1 =\frac{X-\lambda}{\lambda}.$$

Poisson分布では $\operatorname{Var}(X)=\lambda$ だから、

$$\begin{aligned} I_1(\lambda) &=E[U_1(\lambda)^2]\\ &=E\left[\left(\frac{X-\lambda}{\lambda}\right)^2\right]\\ &=\frac{\operatorname{Var}(X)}{\lambda^2} =\frac1\lambda. \end{aligned}$$

独立な $n$ 観測では情報量が足し算されるので、

$$I_1(\lambda)=\frac1\lambda,\qquad I_n(\lambda)=\frac n\lambda.$$

CR下限は $\lambda/n$ です。一方、

$$\operatorname{Var}(\bar X)=\frac{\lambda}{n}.$$

よって $\bar X$ は限界へ達する有効推定量です。

答案で気をつけること： $I_1$ と $I_n$ を混同せず、標本数 $n$ を掛けます。

問題13：正規分布の分散の最尤推定

$X_i\overset{\mathrm{iid}}{\sim}N(0,\theta)$ とし、$\hat\theta=n^{-1}\sum_iX_i^2$ とする。期待値、分散、漸近分布を求めよ。

解答13

方針： $\hat\theta$ は $X_i^2$ の標本平均です。そこで $X_i^2$ の平均と分散を求め、最後に中心極限定理を使います。

$Z_i=X_i/\sqrt\theta\sim N(0,1)$ と置くと、$X_i^2=\theta Z_i^2$ です。$E[Z_i^2]=1$、$E[Z_i^4]=3$ より、

$$E[X_i^2]=\theta,\qquad E[X_i^4]=3\theta^2.$$

したがって、

$$\operatorname{Var}(X_i^2) =3\theta^2-\theta^2=2\theta^2.$$

独立性から、

$$\operatorname{Var}(\hat\theta)=\frac{2\theta^2}{n}.$$

中心極限定理より、

$$\sqrt n(\hat\theta-\theta) \xrightarrow{d}N(0,2\theta^2).$$

答案で気をつけること： $N(0,\theta)$ の第2母数が分散であると確認します。標準偏差なら式が変わります。

問題14：一般の関数に対するCR下限

$T$ が $g(\theta)$ の不偏推定量であるとする。正則条件のもとで分散の下限を示せ。

解答14

方針： 通常のCR不等式で $1$ だった箇所が、$g'(\theta)$ へ変わるだけです。密度の微分を「密度 × スコア」へ置き換えます。

$E[T]=g(\theta)$ を微分すると、

$$g'(\theta) =\frac{\partial}{\partial\theta}E[T] =E[TU_n(\theta)].$$

$E[U_n]=0$ なので $\operatorname{Cov}(T,U_n)=g'(\theta)$ です。Cauchy–Schwarz不等式より、

$$\{g'(\theta)\}^2 \le\operatorname{Var}(T)I_n(\theta).$$

したがって、

$$\operatorname{Var}(T) \ge\frac{\{g'(\theta)\}^2}{I_n(\theta)}.$$

答案で気をつけること： $\theta$ 自身の推定だけなら $g'(\theta)=1$ です。

問題15：最尤推定量の漸近正規性

正則な1母数モデルで、最尤推定量の漸近分布をスコアのTaylor展開から導く道筋を示せ。

解答15

方針： 最尤推定量の誤差を「真値でのスコアの揺れ ÷ 尤度の曲がり具合」と書き直します。その分子にCLT、分母にLLNを使います。

$U_n(\hat\theta)=0$ を真値 $\theta_0$ の周りで展開し、

$$0=U_n(\theta_0)+(\hat\theta-\theta_0)U_n'(\tilde\theta)$$

とします。よって、

$$\sqrt n(\hat\theta-\theta_0) =\frac{U_n(\theta_0)/\sqrt n}{-U_n'(\tilde\theta)/n}.$$

$U_n(\theta_0)=\sum_{i=1}^nU_1^{(i)}(\theta_0)$ は独立な1標本スコアの和で、平均0、分散 $I_1(\theta_0)$ を持ちます。したがって分子へ中心極限定理、分母へ大数の法則を使うと、

$$\frac{U_n(\theta_0)}{\sqrt n} \xrightarrow{d}N\{0,I_1(\theta_0)\}, \qquad -\frac{U_n'(\tilde\theta)}n\xrightarrow{p}I_1(\theta_0).$$

Slutskyの定理より、

$$\sqrt n(\hat\theta-\theta_0) \xrightarrow{d}N\{0,I_1(\theta_0)^{-1}\}.$$

答案で気をつけること： Taylor展開、CLT、LLN、Slutskyのどこを使ったかを分けて書きます。

発展問題16：2母数の情報行列

互いに独立な $X\sim N(\mu,\sigma^2)$ について、1観測の $(\mu,\sigma^2)$ に関するFisher情報行列を求めよ。

解答16

方針： 分散そのものを母数として $v=\sigma^2$ と置き、$\mu,v$ で偏微分します。情報行列の各要素はスコアの積の期待値です。

1観測の対数尤度は

$$\ell(\mu,v) =-\frac12\log(2\pi v)-\frac{(X-\mu)^2}{2v}.$$

スコアは、

$$\begin{aligned} U_\mu&=\frac{\partial\ell}{\partial\mu} =\frac{X-\mu}{v},\\ U_v&=\frac{\partial\ell}{\partial v} =-\frac1{2v}+\frac{(X-\mu)^2}{2v^2}. \end{aligned}$$

$Y=X-\mu$ と置くと、$E[Y]=0$、$E[Y^2]=v$、$E[Y^3]=0$、$E[Y^4]=3v^2$ です。よって対角成分は

$$\begin{aligned} I_{\mu\mu} &=E[U_\mu^2] =\frac{E[Y^2]}{v^2} =\frac1v,\\ I_{vv} &=E\left[\left\{-\frac1{2v}+\frac{Y^2}{2v^2}\right\}^2\right]\\ &=\frac{\operatorname{Var}(Y^2)}{4v^4}\\ &=\frac{3v^2-v^2}{4v^4} =\frac1{2v^2}. \end{aligned}$$

非対角成分は、

$$I_{\mu v} =E[U_\mu U_v] =\frac{E[Y^3]}{2v^3}=0$$

です。$v=\sigma^2$ に戻すと、

$$I(\mu,\sigma^2) = \begin{pmatrix} 1/\sigma^2&0\\ 0&1/(2\sigma^4) \end{pmatrix}.$$

非対角成分が0なので、この母数化では平均と分散は情報の意味で直交しています。

答案で気をつけること： 第2母数を $\sigma$ とするか $\sigma^2$ とするかで情報行列は変わります。

17. 医薬・生命科学の問題

問題1：小標本の薬効反応率

候補化合物を投与した独立な12個体中3個体で反応を認めた。事前分布を $\operatorname{Beta}(2,2)$ とし、反応率の最尤推定値とBayes推定値を求めよ。

解答1

最尤推定値は、

$$\hat p_{\mathrm{ML}}=\frac3{12}=0.25.$$

事後分布は $\operatorname{Beta}(5,11)$ なので、二乗誤差損失でのBayes推定値は、

$$\hat p_B=\frac5{16}=0.3125.$$

少数例なので、事後平均は事前平均0.5の方向へ縮みます。

答案で気をつけること： 「反応した3例」だけでなく総数12も更新に使います。

問題2：希少有害事象率

合計100人年の追跡で4件の有害事象を観測した。発生件数を $\operatorname{Poisson}(100\lambda)$、事前分布を率表示の $\operatorname{Gamma}(1,20)$ とする。事後分布と事後平均を求めよ。

解答2

曝露時間を含む尤度は $\lambda^4e^{-100\lambda}$ に比例します。したがって、

$$\lambda\mid x\sim\operatorname{Gamma}(5,120).$$

事後平均は、

$$E[\lambda\mid x]=\frac5{120}=0.0417$$

件/人年です。

答案で気をつけること： 更新される率母数は観測数ではなく総曝露時間です。単位も添えます。

問題3：2施設の濃度測定を統合する

同じ標準試料の濃度を2施設で測り、$Y_1=9.8$、$Y_2=10.4$、既知標準偏差がそれぞれ1.0、2.0だった。独立・不偏を仮定してBLUEを求めよ。

解答3

まず標準偏差ではなく分散へ直します。

$$\sigma_1^2=1.0^2=1, \qquad \sigma_2^2=2.0^2=4.$$

逆分散は $1$ と $1/4$ です。重みの比は

$$1:\frac14=4:1$$

なので、和が1になるように直すと、

$$w_1=\frac4{4+1}=0.8, \qquad w_2=\frac1{4+1}=0.2.$$

したがって、

$$\begin{aligned} \hat\mu &=0.8\times9.8+0.2\times10.4\\ &=9.92. \end{aligned}$$

逆分散の式へ直接代入しても同じです。

$$\hat\mu =\frac{1\cdot9.8+(1/4)\cdot10.4}{1+1/4} =9.92.$$

精度の高い施設1へ強く重み付けされます。

答案で気をつけること： 標準偏差の逆数ではなく、分散の逆数を重みにします。また、同じ真の濃度を測っているという前提が崩れるなら、BLUEだけで施設差を処理してはいけません。

問題4：濃度反応曲線の最尤推定

4母数logisticモデルでEC50を推定するとき、最尤推定値だけを報告する危険性を3つ挙げ、追加すべき評価を答えよ。

解答4

危険性は、例えば次の3つです。

1. EC50とHill係数が強く相関し、尤度が平らな場合がある
2. 上限・下限が観測範囲外だとEC50が不安定になる
3. 濃度によって分散が変わると、等分散正規尤度が不適切になる

プロファイル尤度、bootstrap、残差図、母数間相関、独立実験間の再現性を併記します。

答案で気をつけること： ウェル数を独立実験数として扱わず、プレートや実験日の階層を確認します。

問題5：個体内に細胞があるデータ

対照3個体、処置3個体から各個体1000細胞を測定した。6000細胞を独立として平均差を推定する問題点と、適切なモデル方針を説明せよ。

解答5

同じ個体の細胞は共通の生物学的背景を持つため独立ではありません。6000を標本数にすると標準誤差を過小評価する擬似反復になります。

個体を実験単位とし、個体ごとの集約値を解析するか、細胞を個体内にネストした混合モデル・階層Bayesモデルを使います。

答案で気をつけること： 観測単位の細胞数と、独立な実験単位の個体数を分けて書きます。

問題6：集団薬物動態と経験Bayes

母集団解析で得られる個体別経験Bayes推定値が、観測の少ない個体ほど母集団平均へ近づく理由を説明せよ。

解答6

個体データが少ないと個体尤度の情報量が小さく、母集団分布という事前情報の相対的重みが大きくなります。そのため個体推定値は母集団平均へ強く縮みます。

観測が多い個体では尤度の情報量が増え、個体データ側の重みが大きくなります。

答案で気をつけること： 縮小を「補正」とだけ呼ばず、事前情報と個体尤度の情報量の釣り合いとして説明します。

問題7：RNA-seqの分散推定

遺伝子ごとに少数標本しかないとき、各遺伝子の分散を完全に別々に推定するより、経験Bayesで情報共有する利点と注意点を述べよ。

解答7

多数遺伝子から分散の全体傾向を学び、極端に不安定な遺伝子別分散を縮小できるため、MSEと検定の安定性が改善します。

一方、全遺伝子が似た分散構造を持つという仮定、超母数推定の不確実性、強い外れ遺伝子への頑健性を確認する必要があります。

答案で気をつけること： 細胞数やread数が多くても、生物学的replicateが少ない問題は解消されません。

問題8：L2正則化とMAP

正規誤差の線形回帰で係数 $\beta$ に平均0の正規事前分布を置くと、MAP推定がL2正則化に対応することを説明せよ。

解答8

正規尤度の負の対数は残差平方和に比例し、正規事前分布の負の対数は $\|\beta\|_2^2$ に比例します。したがって事後最大化は、

$$\operatorname*{minimize}_\beta \left\{ \sum_i(y_i-x_i^\mathsf T\beta)^2 +\lambda\|\beta\|_2^2 \right\}$$

と同値です。

答案で気をつけること： MAPの点推定と、事後分布全体を用いたBayes予測を区別します。

問題9：多バイオマーカーの情報行列

2つのバイオマーカー効果 $\theta_1,\theta_2$ の情報行列の非対角成分が大きいとき、推定上何が起こるか説明せよ。

解答9

2つの母数の尤度方向が強く結び付き、片方を変えた影響をもう片方で補えるため、個々の効果を分離しにくくなります。情報行列の逆の対角成分が大きくなり、標準誤差が増えることがあります。

追加デザイン、直交化、事前情報、複合指標の利用を検討します。

答案で気をつけること： 非対角成分そのものを相関係数と断定せず、逆行列から推定量の共分散を評価します。

問題10：毒性予測の不確実性

毒性分類モデルの出力確率を「真の毒性確率」と解釈する前に、推定の観点から確認すべき項目を挙げよ。

解答10

独立な外部検証での校正、クラス不均衡、学習分布と適用先の差、モデル・母数・データ由来の不確実性、閾値ごとの損失を確認します。

深層学習のsoftmax値は自動的に校正された事後確率にはなりません。温度スケーリング、bootstrap、ensembleなどで不確実性と校正を評価します。

答案で気をつけること： 識別性能のAUROCだけでは、確率推定の正しさや臨床的損失は評価できません。

18. 章のまとめ

• 十分統計量は、母数の情報を失わずに標本を圧縮する
• 因子分解定理は十分性を尤度の分解から判定する
• 指数型分布族では自然母数、十分統計量、対数分配関数がつながる
• モーメント法は母集団モーメントと標本モーメントを等置する
• 最尤法は観測データに対する尤度を最大化する
• Bayes法は事前分布と尤度を事後分布へ統合する
• 推定量は不偏性だけでなく、分散、MSE、一致性で比べる
• Fisher情報量は尤度の曲率と推定精度を表す
• Cramér–Rao不等式は正則な不偏推定の分散下限を与える
• 最尤推定量は正則条件下で一致・漸近正規・漸近有効になる
• 階層構造のある医薬・生命科学データでは、実験単位を守って部分プーリングする