ダメ出し：ifelse は遅い

Rでフィボナッチ数！ http://d.hatena.ne.jp/mickey24/20080914/1221334119
syou6162先生が一行でやってくれました http://d.hatena.ne.jp/mickey24/20080914/1221334177
Rでフィボナッチ数！(動的計画法版) http://d.hatena.ne.jp/mickey24/20080914/1221336203

再帰関数として書いたもの　これは，おそい
fib1 <- function(n) {
  if (n == 0) { return(0) }
  if (n == 1) { return(1) }
  fib1(n-1) + fib1(n-2)
}

一行野郎で書いたもの　これは，作者は心外だろうが，大変大変遅い
なぜおそいかというと，それは，ifelse を使っているから
ifelse は，条件によらず二通りの計算をして，条件によりどちらかを選ぶ　だから，とってもとっても遅い
fib2 <- function(n) {
    ifelse (n > 1, Recall(n-1) + Recall(n-2), ifelse(n == 1, 1, 0))
}

同じ一行にするのでも，if else を使えば func1 と基本的に同じなのだから，ifelse ほど遅くはない
fib3 <- function(n) {
    if (n > 1) Recall(n-1) + Recall(n-2) else if (n == 1) 1 else 0
}

動的計画法版のプログラムだが，ベクトルを使う必要は余り感じられない
fib4 <- function(n) {
  if (n == 0) { return(0) }
  x <- numeric(max(n,2))
  x[1] <- x[2] <- 1
  m <- 3
  while (m <= n) {
    x[m] <- x[m-1] + x[m-2]
    m <- m + 1
  }
  x[n]
}

普通の変数を使って書けば，最速!!
fib5 <- function(n) {
    if (n == 0) return(0)
    if (n == 1) return(1)
    a <- 0
    b <- 1
    for (i in 2:n) {
        c <- a + b
        if (i == n) return(c)
        a <- b
        b <- c
    }
}

> n <- 15
> benchmark(fib1(n), fib2(n), fib3(n), fib4(n), fib5(n),
+           columns=c("test", "elapsed", "relative", "replications"),
+           replications=200)
     test elapsed    relative replications
1 fib1(n)   0.751  125.166667          200
2 fib2(n)   6.229 1038.166667          200
3 fib3(n)   0.803  133.833333          200
4 fib4(n)   0.014    2.333333          200
5 fib5(n)   0.006    1.000000          200

ダメ出し：重回帰分析の標準誤差の推定値(Residual standard error)

C.P. ロバート， G. カセーラ 著
石田基広，石田和枝
R によるモンテカルロ法入門

23ページに，「lm() 関数は（中略）意外なことに標準誤差の推定値がリスト要素として含まれておらず（中略）このすうちそのものは次のようにして求める必要があります」

これは，summary.lm 関数が返す $sigma である。わざわざ計算しなくても良い。
回帰の分散分析の F 統計量と自由度も lm 関数の戻り値にはないが，summary.lm 関数の戻り値には含まれる。しかし，その P 値は自分で計算しないといけない。

> x <- 1:10
> y <- c(9, 7, 1, 4, 5, 10, 6, 2, 8, 3)
> a <- lm(y ~ x)
> (b <- summary(a))

Call:
lm(formula = y ~ x)

Residuals:
    Min      1Q  Median      3Q     Max
-4.9697 -1.7500  0.0939  2.2015  4.5939

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept)   6.5333     2.1547   3.032   0.0163 *
x            -0.1879     0.3473  -0.541   0.6032  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.154 on 8 degrees of freedom
Multiple R-squared: 0.0353,    Adjusted R-squared: -0.08529
F-statistic: 0.2927 on 1 and 8 DF,  p-value: 0.6032

> (u <- sqrt(sum(a$res^2) / a$df))
[1] 3.154122
> (v <- b$sigma) # これですよ!
[1] 3.154122
> all(u == v)
[1] TRUE

> b$fstatistic
    value     numdf     dendf
0.2927201 1.0000000 8.0000000

> (p <- pf(b$fstatistic[1], b$fstatistic[2],  b$fstatistic[3], lower.tail=FALSE))
    value
0.6032176

barplot 関数の戻り値

barplot 関数は，描いた棒の中心 x 座標を返す。

座標軸，目盛り等を描くときにはそれを用いればよい。

> x <- dbinom(0:6, 6, 0.1)
> a <- barplot(x)
> str(a)
 num [1:7, 1] 0.7 1.9 3.1 4.3 5.5 6.7 7.9
> axis(1, a, 0:6)

ダメ出し：小細工せずに素直に（あるいは，効果を確かめてから） その2

驚くべきことに（??）

    w <- matrix(a, 2)[2,]

は

    v <- a[1:(n%/%2)*2L]                        # 最速！！

と同じくらいに速い。