読者です 読者をやめる 読者になる 読者になる

餡子付゛録゛

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

「Rと手作業で覚える最尤法」のコードに関して補足

どうもコメント欄で粗雑なコードで御迷惑をおかけしましたと書いたら、『裏 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マニアっぽいので、エレガントなコードを書きたい人はそちらを読まれる事をお勧めします。