どうもコメント欄で粗雑なコードで御迷惑をおかけしましたと書いたら、『裏 RjpWiki』と言うブログで「Rと手作業で覚える最尤法」に関して、『ある意味「えれがんと」だなあ』と皮肉を言われています（ダメ出し：簡明直截なプログラムを書こう）。
ダメ出しを受けているわけですが、どうも元のコードがこれまたエレガントで無かった為か、速度のためか、“直截簡明”にして頂いたのにループ内で行列式の作成が大量に残ってしまっています。この二つをループ内で構成するのは、見通しの面で良くないかも知れません。ちょっと整理をしてみました。
# サンプルを初期化
x <- c(11, 13, 23)
n <- length(x)
# 対数尤度関数のグラディエント・ベクトル
f <- expression(matrix(c(
-sum(x-m[1])/m[2],
-n/(2*m[2]) + sum( (x-m[1])^2)/(2*m[2]^2 )
), 2, 1))
# muとs2でヤコビアンを作っておく
g <- expression(matrix(c(
n/m[2],
sum(m[1] - x)/m[2]^2,
sum(x - m[1])/m[2]^2,
n/(2*m[2]^2) - sum( (x - m[1])^2 )/m[2]^3
), 2, 2))
# 平均値が10、分散10で初期化
m <- matrix(c(10, 10), 2, 1)
# ループする
for(i in 1:10){
m <- m - solve(eval(g), eval(f))
print(sprintf("[%d] (mu,s2)=(%f,%f)", i, m[1], m[2]))
}
# 計算結果を表示
print(sprintf("平均%2.3f、分散%2.3f（標準偏差%2.3f）の正規分布", m[1], m[2], sqrt(m[2])))
速度的には修正して頂いたコードよりも2割程度遅くなるわけですが、ループ処理内で行列式の作成に注意を払わなくて済み、さらに説明手順を保持できるので、直観性は維持できているのではないでしょうか。
このブログのコードは大雑把な説明用にできているので、「裏 RjpWiki」の中の人を混乱させてしまっているのかも知れません。「裏 RjpWiki」の中の方はRマニアっぽいので、エレガントなコードを書きたい人はそちらを読まれる事をお勧めします。
