Rと手作業で覚える最尤法

OLSより進んだ統計手法で最初に覚えるのは最尤法だと思います。大半の人はツールとして知っていて、あまり中身を意識していない気がするのですが、「尤度」の説明無しで『尤度が最大になるパラメーターを求める方法』と言う説明が横行しているのは、問題があるかも知れません。
最尤法は、ある分布から観測値が取り出されたとして、“そうなる確率”が最も高くなるように分布の具体的な形状を決めるやり方です。“そうなる確率”を尤度と言います。こう書くと易しい事なのか難しい事なのか判別もつかないと思うので、実際に最尤法を解いてみましょう。
まず、何も考えずにトライ＆エラーで最尤法を試みるやり方を説明した後に、教科書的な最尤法の解法を説明します。

1. 何も考えずにトライ＆エラーで最尤法を試みる

ある正規分布から値を3つ取り出したら、11 13 23だったとしましょう。このサンプルが“もっともらしい”正規分布の平均と分散は何になるでしょうか。“もっともらしい”程度を各観測値が抽出される確率を乗じた同時確率だと定義し、それを尤度と呼ぶことにします。
Rで平均15、分散40の正規分布の尤度を計算すると、次のようになります。

# 観測値のセット
x <- c(11, 13, 23)
# 平均15、分散40の正規分布から観測値それぞれが観測される確率を計算
p <- dnorm(x, mean=15, sd=sqrt(40))
# 確率を乗じて、同時確率にする
print(paste("尤度:", prod(p)))

平均と分散を山勘で決めて尤度を求めていきましょう。

平均	分散	11の確率	13の確率	23の確率	尤度（＝同時確率）
15	40	0.053	0.061	0.027	0.000085708
15	41	0.053	0.060	0.027	0.000084767
15	42	0.052	0.059	0.027	0.000083807
16	40	0.048	0.058	0.032	0.000088955
16	41	0.048	0.057	0.032	0.000087899
16	42	0.047	0.056	0.032	0.000086830
17	40	0.042	0.053	0.038	0.000085708
17	41	0.042	0.053	0.038	0.000084767
17	42	0.042	0.052	0.038	0.000083807

ここでの尤度は(11の確率)×(13の確率)×(23の確率)になる事に注意してください。平均が15で分散が40の正規分布であれば、11は0.053、13は0.061、23は0.027の確率で発生するので、尤度はそれらを乗じて0.000085708になります。
平均が16、分散が40の所で尤度が最大になりました。試した範囲では、これが“もっともらしい”正規分布の平均と分散となります。観測されたサンプルの値は固定（11、13、23）になるので、尤度関数は平均と分散の関数になります。

2. 尤度関数、対数尤度関数、一階条件から最尤法を試みる

実際は手当たり次第にパラメーター（平均と分散）の組み合わせを試すわけにはいかないので、尤度関数の対数をとって対数尤度関数にした後に、対数尤度関数をパラメーターでそれぞれ微分して連立方程式にした後に、数値解析アルゴリズムを使ってせっせと計算します。
尤度関数は平均と分散をパラメーターとする同時確率となりますから、以下のようになります（「経済統計補足資料#27：最尤推定量の導出」から転載）。要するに正規分布による同時確率の計算式ですね。

$\mathfrak{L}(\mu, \sigma^2) = \frac{1}{\sqrt{2 \pi \sigma^2}}\exp \big( - \frac{(x_1 - \mu)^2}{2\sigma^2} \big ) \cdot \frac{1}{\sqrt{2 \pi \sigma^2}}\exp \big( -\frac{(x_2 - \mu)^2}{2\sigma^2} \big) \cdots \frac{1}{\sqrt{2 \pi \sigma^2}}\exp \big( -\frac{(x_n - \mu)^2}{2\sigma^2} \big) \\ = \big( \frac{1}{\sqrt{2\pi}} \big)^n \big( \frac{1}{\sqrt{\sigma^2}} \big)^n \exp{ \big( - \frac{1}{2\sigma^2} \sum^n_{i=1} (x_i-\mu)^2 \big) }$

このままだと計算が難しいので対数化します。

$\log{ \mathfrak{L}(\mu, \sigma^2)} = \log{ (2\pi)^{-\frac{n}{2}}} + \log{ (\sigma^2)^{-\frac{n}{2}}} + \log{ \exp \big( -\frac{1}{2\sigma^2} \sum^n_{i=1} (x_i - \mu)^2 \big) } \\ = - \frac{n}{2} \log{2\pi} - \frac{n}{2} \log{\sigma^2} - \frac{1}{2\sigma^2} \sum^n_{i=1} (x_i - \mu)^2$

対数化したままだとコンピューターは計算できないので、偏微分して一階の条件を整理します。

$\frac{ \partial \log{ \mathfrak{L}(\mu, \sigma^2)} }{ \partial \mu } = (-1) \cdot 2 \cdot \frac{1}{2\sigma^2} \sum^n_{i=1} (x_i - \mu) = - \frac{1}{\sigma^2} \sum^n_{i=1} (x_i - \mu)$

$\frac{ \partial \log{ \mathfrak{L}(\mu, \sigma^2)} }{ \partial \sigma^2 } = \frac{n}{2\sigma^2} + \frac{1}{2 ( \sigma^2 )^2 } \sum^n_{i=1} (x_i - \mu)^2$

この連立方程式をニュートン・ラフソン法のような数値解析アルゴリズムで解きます。正規分布の場合は、整理すると手計算可能なぐらいシンプルな方程式になりますが、試しにコードを書いてみましょう。

# サンプル
x <- c(11, 13, 23)
# x <- rnorm(300, mean=15, sd=10) 等でも計算できる
n <- length(x)
# 連立方程式を設定
f1 <- expression(-sum(x-mu)/s2)
f2 <- expression(-n/(2*s2) + sum((x-mu)^2)/(2*(s2^2)))
# muとs2で一階微分を作っておく
g11 <- expression(n/s2)
g12 <- expression(sum(x - mu)/(s2^2))
g21 <- expression(sum(mu - x)/(s2^2))
g22 <- expression(n/(2*(s2^2)) - sum( (x - mu)^2 )/(s2^3))

# 初期値を設定
mu <- 10 # 正規分布の平均値
s2 <- 10 # 正規分布の分散
for(i in 1:10){
	# (mu, s2)を2x1行列にする
	m <- matrix(c(mu, s2), 2, 1)
	# (mu, s2)で評価したf1とf2の値を2x1行列にする
	f <- matrix(c(eval(f1), eval(f2)), 2, 1)
	# ヤコビアンに(mu, s2)を代入した行列を作る
	j <- matrix(c(eval(g11), eval(g21), eval(g12), eval(g22)), 2, 2)
	# 行列mから、ヤコビアンの逆行列と評価式を乗じたものを引く
	m <- m - solve(j)%*%f
	# 行列mを(x, y)に展開しておく
	mu <- m[1]; s2 <- m[2]
	print(sprintf("[%d] (mu,s2)=(%f,%f)", i, mu, s2))
}
print(sprintf("平均%2.3f、分散%2.3f（標準偏差%2.3f）の正規分布", mu, s2, sqrt(s2)))

実行すると以下の表のように9順目ぐらいで収束します。初期値を例えば0付近にすると収束しなくなるので注意してください。

ステップ	平均	分散	11の確率	13の確率	23の確率	尤度（＝同時確率）
1	17.5	6.8	0.007	0.034	0.016	0.000003844
2	16.3	9.2	0.028	0.072	0.012	0.000023852
3	15.9	12.8	0.043	0.080	0.016	0.000054561
4	15.8	17.3	0.050	0.077	0.021	0.000080688
5	15.7	21.9	0.052	0.072	0.025	0.000093919
6	15.7	25.7	0.051	0.069	0.028	0.000097562
7	15.7	27.3	0.051	0.067	0.029	0.000097937
8	15.7	27.6	0.051	0.067	0.029	0.000097943
9	15.7	27.6	0.051	0.067	0.029	0.000097943

対数尤度関数の導出は統計学の論文一本分の学術的貢献になるので、大半の応用経済学者は単に対数尤度関数をコピーして利用しています。Rの場合は対数尤度関数を導出したら、後はnlm()関数がやってくれます。さらにプロビット/ロジット・モデルなどの利用例が多いモデルは、パッケージが用意されているので最尤法も対数尤度関数も意識する必要はありません。

3. 制約条件をつけて多変量回帰分析にする場合

上の例では制約条件がありませんでした。つまり重回帰分析の係数は無く、単に平均と分散を求めただけになります。重回帰分析の場合でも、対数尤度関数のx-muの部分がy-b0-b1*x1-b2*x2などと置き換わるだけなので、概念的には簡単です。ただし推定するパラメーター数と同じだけ一階条件の式が増え、ヤコビアンはその二乗の大きさになる事に注意してください。
とてもでは無いですが、上の例のようにパッケージ無しで計算するのには無理があります。

4. t値、P値などの計算

テキストを見れば書いてありますが、最尤法のヘッセ行列の逆行列の対角成分は推定したパラメーターの分散に等しくなります。ここで言うヘッセ行列は、最尤法の対数尤度関数の二階微分を推定量で評価したものです。

SEs <- sqrt(abs(diag(solve(j))))

それぞれの係数を、それぞれの標準誤差で割ると、帰無仮説（H₀）が0のt値が出て来ます。

tv <- m / SEs

両側検定を行えば、パッケージが表示するP値になります。自由度は(サンプルサイズ−パラメーターの数)になります。

tp <- 2*(1 - pt(abs(tv), n - length(tv)))
sprintf("t値%2.3f P値%2.3f",tv,tp)

平均値は12.2%有意、分散は34.6%有意です。観測数が3しかないので、統計的有意性が高くないと言う事のようです。2変数ですが、10から20のサンプルは必要なようです。

5. 終わりに

尤度関数や対数尤度関数を間違えると正しくない推定になったり、はたまたその微分を間違えると収束しないと言う現象が発生します。実用上は、なるべく統計解析パッケージを利用した方が良いでしょう。しかし裏側で何をやっているかを知っておく事は、大切な事だと思います。

餡子付゛録゛

ソフトウェア開発ツールの便利な使い方を紹介。