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
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 関数は,描いた棒の中心 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)
驚くべきことに(??)
w <- matrix(a, 2)[2,]
は
v <- a[1:(n%/%2)*2L] # 最速!!
と同じくらいに速い。