ダメ出し：簡明直截なプログラムを書こう

このブログは，プログラムの速度を追求するのではなく，「えれがんと」なプログラミングを目指しているので誤解ないように。「えれがんと」が何を意味するかは範囲は広い。

で，「2012-01-10  Rと手作業で覚える最尤法」の「2. 尤度関数、対数尤度関数、一階条件から最尤法を試みる」について

このプログラムでは，導関数の定義を expression で行っておいてから，使うときに eval で評価している。なるほど，このようにするとある意味「えれがんと」だなあ。

でも，直接指定しても「えれがんと」でなくなるとも見えない。結果としてわずかだが速いし。

なお，連立方程式を解くとき，solve で逆用列を求めて定数項と行列掛け算するのではなく，定数項を solve の引数に渡す方がよいというのは，数値演算の定石。

もとのプログラム（整形済み）

> system.time({for (i in 1:10000) {
+ x <- c(11, 13, 23)
+ n <- length(x)
+ f1 <- expression(-sum(x - mu) / s2)
+ f2 <- expression(-n / (2 * s2) + sum((x - mu) ^ 2) / (2 * (s2 ^ 2)))
+ 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) {
+     m <- matrix(c(mu, s2), 2, 1)
+     f <- matrix(c(eval(f1), eval(f2)), 2, 1)
+     j <- matrix(c(eval(g11), eval(g21), eval(g12), eval(g22)), 2, 2)
+     m <- m - solve(j) %*% f
+     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)))
+ }})
   ユーザ   システム       経過  
    14.484      0.062     14.270

書き換えたプログラム

> system.time({for (i in 1:10000) {
+ x <- c(11, 13, 23)
+ n <- length(x)
+ mu <- 10
+ s2 <- 10
+ for (i in 1:10) {
+     m <- c(mu, s2)
+     f <- c(-sum(x - mu)/s2,
+            -n/(2 * s2) + sum((x - mu)^2)/(2 * (s2^2)))
+     j <- matrix(c(n/s2,
+                   sum(mu - x)/(s2^2),
+                   sum(x - mu)/(s2^2),
+                   n/(2 * (s2^2)) - sum((x - mu)^2)/(s2^3)), 2, 2)
+     m <- m - solve(j, f)
+     mu <- m[1]
+     s2 <- m[2]
+ #    cat(sprintf("[%d] (mu, s2) = (%f, %f)\n", i, mu, s2))
+ }
+ #cat(sprintf("平均%.3f、分散%.3f（標準偏差%.3f）の正規分布\n",
+ #    mu, s2, sqrt(s2)))
+ }})
   ユーザ   システム       経過  
     8.878      0.022      8.879