2011-11-06 ２次計画問題
http://d.hatena.ne.jp/ryamada22/20111106/1320544702
segments は，ベクトル化されている
# for(i in 1:N){
# segments(i,D[i,1],i,D[i,2],col=3)
#}
ではなく
segments(1:N, D[, 1], 1:N, D[, 2], col=3)
しかし，
for(i in 1:N){
for(j in 2:M){
segments(i,D[i,1],i,D[i,j],col=1+j)
}
}
を
segments(1:N,D[,1],1:N,D[,2:M],col=1+2:M)
だけにしてしまうと，segments の色が違ってしまうなあ。
R 2.13のコンパイラをつかってみる(3)
http://ito-hi.blog.so-net.ne.jp/2011-07-09-2
> set.seed(1234)
> n <- 10000000 # 速すぎて速度が測れないので多めに設定
>
> x <- rnorm(n, 0, 1)
>
> # ctest2 は test2 をコンパイルしたもの
> test2 <- function(x) {
+ d <- vector("numeric", length(x) - 1)
+ for (i in 1:(length(x) - 1)) {
+ d[i] <- exp(x[i + 1] - x[i])
+ }
+ mean(d)
+ }
>
> ctest2 <- cmpfun(test2)
>
> system.time(a <- ctest2(x)) # for を使ったコンパイラ版（4 つの中では最速）
ユーザ システム 経過
11.379 0.151 11.585
> a
[1] 2.717731
>
> # test2 をベクトル演算を使って高速化
> test3 <- function(x) {
+ n <- length(x)
+ mean(exp(x[-1]-x[-n]))
+ }
> system.time(ans3 <- test3(x))
ユーザ システム 経過
0.669 0.268 0.930 # これでも，17 倍速
> ans3
[1] 2.717731
>
> # test2 を インラインにする
> library(inline)
> src <- '
+ Rcpp::NumericVector x(X);
+ int n = x.length();
+ double mean = 0;
+ for (int i = 0; i < n-1; i++) {
+ mean +=exp(x(i+1)-x(i));
+ }
+ return wrap(mean/(n-1));
+ '
> test4 <- cxxfunction(signature(X="numeric"), src, plugin="Rcpp")
> system.time(ans4 <- test4(x))
ユーザ システム 経過
0.251 0.001 0.250 # 45 倍速
> ans4
[1] 2.717731
R 2.13.0のコンパイラをつかってみる (2) Cλ
http://ito-hi.blog.so-net.ne.jp/2011-04-21
Morisita の Cλ を求めるのに，行ごとの組合せで n.cmm × n.cmm 回 clambda 関数を呼んでいるが，無駄。
対称行列なのだから，せめて半分だけ計算する。
以下，原著者によるプログラムのパフォーマンス
> system.time(a1 <- mclambda(cmm))
ユーザ システム 経過
10.360 0.136 10.423
> system.time(a2 <- mclambda2(cmm)) # コンパイル版
ユーザ システム 経過
7.993 0.188 8.127
> system.time(a3 <- c.mclambda(cmm)) # コンパイル版
ユーザ システム 経過
7.093 0.278 7.333
ベクトル演算を使えば，コンパイル無しで，コンパイルしたものより 100 倍速くなる。
以下だけで，480×480 の組合せの Cλ が一挙に得られる。
> system.time({
+ n1n2 <- rowSums(cmm)
+ l1l2 <- rowSums(cmm*(cmm-1))/(n1n2*(n1n2-1))
+ cc <- cmm%*%t(cmm)
+ l1pl2 <- outer(l1l2, l1l2, "+")
+ n1mn2 <- outer(n1n2, n1n2, "*")
+ ans <- 2*cc/(l1pl2*n1mn2)
+ })
ユーザ システム 経過
0.072 0.011 0.083 # 100 倍速!!
当然ながら，結果は等しい
> all.equal(a1, ans)
[1] TRUE
> all.equal(a2, ans)
[1] TRUE
> all.equal(a3, ans)
[1] TRUE
移動平均は stats::filter で行える。 これを，inline で書いて見ようというのが今回のお題。
NA は NA_REAL を使う（IntegerVector なら，NA_INTEGER。他に，Inf, -Inf は R_PosInf, R_NegInf。）
ポインターは double *pz; ... pz = REAL(z); のように使う。（ので，いいのかな?）
かなりいい加減に書いたプログラムだけど，filter より 6 倍ほど速い。もっとも，要素数が 1000000 と大きいのに，実行時間が 0.15 秒というのだから，ほとんど意味がない。
src <- '
Rcpp::NumericVector x(X);
int n = x.length();
Rcpp::NumericVector z(n, NA_REAL);
double *pz;
int m = as<int>(M);
int i, j;
pz = REAL(z)+m/2;
for (i = 0; i <= n-m; i++) {
double ans = 0;
for (j = 0; j < m; j++) {
ans += x[i+j];
}
*pz++ = ans/m;
}
return wrap(z);
'
moving.average <- cxxfunction(signature(X="numeric", M="integer"), src, plugin="Rcpp")
実行例
> set.seed(777)
> x <- rnorm(1000000)
> system.time(a <- filter(x, rep(1/7, 7)))
ユーザ システム 経過
0.150 0.011 0.159
> system.time(b <- moving.average(x, 7))
ユーザ システム 経過
0.026 0.000 0.027
> all.equal(as.vector(a), b)
[1] TRUE
