64 bit integer 登場

int64 パッケージ

> library(int64)
> a <- as.int64(1)
> for (i in 1:300) {a <- a*i; print(a)}
[1] 1
[1] 2
[1] 6
[1] 24
[1] 120
[1] 720
[1] 5040
[1] 40320
[1] 362880
[1] 3628800
[1] 39916800
[1] 479001600
[1] 6227020800
[1] 87178291200
[1] 1307674368000
[1] 20922789888000
[1] 355687428096000
[1] 6402373705728000
[1] 121645100408832000
[1] 2432902008176640000
[1] NA 　　オーバーフロー

gmp パッケージ

> library(help=gmp)
> a <- as.bigz(1)
> for (i in 1:21) { a <- a*i; print(a) }
[1] "1"
[1] "2"
[1] "6"
[1] "24"
[1] "120"
[1] "720"
[1] "5040"
[1] "40320"
[1] "362880"
[1] "3628800"
[1] "39916800"
[1] "479001600"
[1] "6227020800"
[1] "87178291200"
[1] "1307674368000"
[1] "20922789888000"
[1] "355687428096000"
[1] "6402373705728000"
[1] "121645100408832000"
[1] "2432902008176640000"
[1] "51090942171709440000"
[1] "1124000727777607680000"
[1] "25852016738884976640000"

単なる numeric

> options(digits=22)
> a <- 1
> for (i in 1:23) { a <- a*i; print(a) }
[1] 1
[1] 2
[1] 6
[1] 24
[1] 120
[1] 720
[1] 5040
[1] 40320
[1] 362880
[1] 3628800
[1] 39916800
[1] 479001600
[1] 6227020800
[1] 87178291200
[1] 1307674368000
[1] 20922789888000
[1] 355687428096000
[1] 6402373705728000
[1] 121645100408832000
[1] 2432902008176640000
[1] 51090942171709440000
[1] 1124000727777607680000
[1] 25852016738884978212864 　　オーバーフロー！！

等式を成り立たせよう

等式を成り立たせよう
さっきのと同じように解ける。ユニークな解は1つだけど。

> library(e1071)
> invisible(apply(permutations(9), 1, function(i) {
+     st <- paste(c("", "/", "", "+", "/", "", "+", "/", ""), i, collapse="", sep="")
+     if (eval(parse(text=st)) == 1) cat(st, "\n")
+ }))
5/34+7/68+9/12
7/68+5/34+9/12
7/68+9/12+5/34
5/34+9/12+7/68
9/12+5/34+7/68
9/12+7/68+5/34

別解

> library(e1071)
> invisible(apply(permutations(9), 1, function(i) {
+     d1 <- i[2]*10+i[3]
+     d2 <- i[5]*10+i[6]
+     d3 <- i[8]*10+i[9]
+     if ((i[1]*d2*d3+i[4]*d1*d3+i[7]*d1*d2)== d1*d2*d3) cat(i, "\n")
+ }))
5 3 4 7 6 8 9 1 2
7 6 8 5 3 4 9 1 2
7 6 8 9 1 2 5 3 4
5 3 4 9 1 2 7 6 8
9 1 2 5 3 4 7 6 8
9 1 2 7 6 8 5 3 4

式を成り立たせる問題

式を成り立たせる問題。簡単すぎる。×を*，÷を/，=を==とする。

> library(e1071)
> op <- c("+", "-", "*", "/", "==")
> invisible(apply(permutations(5), 1, function(i) {
+     st <- paste(paste(5:1, op[i], collapse="", sep=""), "0", sep="")
+     if (eval(parse(text=st))) cat(st, "\n")
+ }))
5-4*3/2+1==0
5==4*3/2-1+0

pbirthday の改訂

R 2.14.0 から，
pbirthday() and qbirthday() now use exact calculations for coincident = 2.
になった。

R 2.13.1 では，
pbirthday() and qbirthday() did not implement the algorithm exactly as given in their reference and so were unnecessarily inaccurate.
などとうそぶいており，

ずっとまえから，近似解だよといいながら不適切な関数を提供し続けてきた。やっと悔い改めたと言うことか。頑固なんだから。

誕生日のパラドックスは，R 1.9.0 あたりかな。

既約分数クイズ

既約分数クイズ にあるもの。

0/q, q/q は自明なので，除く。いかにもRらしい解??。ただし，メモリ節約も最適化も何にも考えない。既存関数を利用し，ひたすら短く難解なプログラムを目指す。

> func <- function(d0)
+ {
+     euclid <- function(m, n)
+     {
+         m0 <- m
+         n0 <- n
+         while ((temp <- n %% m) != 0) {
+             n <- m
+             m <- temp
+         }
+         return(c(m0/m, n0/m))
+     }
+     a <- unique(data.frame(t(combn(d0, 2, function(x) euclid(x[1], x[2])))))
+     apply(a[order(a[,1]/a[,2]), 1:2], 1, function(x) cat(sprintf("%d/%d  ", x[1], x[2])))
+     cat("\n")
+ }
> func(9)
1/9  1/8  1/7  1/6  1/5  2/9  1/4  2/7  1/3  3/8  2/5  3/7  4/9  1/2  5/9  4/7  3/5  5/8  2/3 
5/7  3/4  7/9  4/5  5/6  6/7  7/8  8/9 

別解（さらに横着）

> library(gmp)
> func <- function(d0)
+ {
+     euclid <- function(m, n)
+     {
+         k <- as.numeric(gcd.bigz(m, n))
+         return(c(m/k, n/k))
+     }
+     a <- unique(data.frame(t(combn(d0, 2, function(x) euclid(x[1], x[2])))))
+     apply(a[order(a[,1]/a[,2]), 1:2], 1, function(x) cat(sprintf("%d/%d  ", x[1], x[2])))
+     cat("\n")
+ }
> func(9)
1/9  1/8  1/7  1/6  1/5  2/9  1/4  2/7  1/3  3/8  2/5  3/7  4/9  1/2  5/9  4/7  3/5  5/8  2/3
5/7  3/4  7/9  4/5  5/6  6/7  7/8  8/9 

さらに短く複雑に

> library(gmp)
> func <- function(d0)
+ {
+     a <- unique(data.frame(t(combn(d0, 2, function(x) {k <- as.numeric(gcd.bigz(x[1], x[2]));
+          return(c(x[1]/k, x[2]/k))}))))
+     apply(a[order(a[,1]/a[,2]), 1:2], 1, function(x) cat(sprintf("%d/%d  ", x[1], x[2])))
+     cat("\n")
+ }
> func(9)
1/9  1/8  1/7  1/6  1/5  2/9  1/4  2/7  1/3  3/8  2/5  3/7  4/9  1/2  5/9  4/7  3/5  5/8  2/3
5/7  3/4  7/9  4/5  5/6  6/7  7/8  8/9 

sd(データフレーム) は廃止-その2

なぜ廃止されたのかな？ストイックな理由なんだろうか。

現在の sd は

> sd
function (x, na.rm = FALSE)
{
    if (is.matrix(x)) {
        msg <- "sd(<matrix>) is deprecated.\n Use apply(*, 2, sd) instead."
        warning(paste(msg, collapse = ""), call. = FALSE, domain = NA)
        apply(x, 2, sd, na.rm = na.rm)
    }
    else if (is.vector(x))
        sqrt(var(x, na.rm = na.rm))
    else if (is.data.frame(x)) {
        msg <- "sd(<data.frame>) is deprecated.\n Use sapply(*, sd) instead."
        warning(paste(msg, collapse = ""), call. = FALSE, domain = NA)
        sapply(x, sd, na.rm = na.rm)
    }
    else sqrt(var(as.vector(x), na.rm = na.rm))
}
<bytecode: 0x1083677b8>
<environment: namespace:stats>

なんだから，warning だけを除いて，従来通りにすればよいのではないか？

sd
function (x, na.rm = FALSE)
{
    if (is.matrix(x))
        apply(x, 2, sd, na.rm = na.rm)
    else if (is.vector(x))
        sqrt(var(x, na.rm = na.rm))
    else if (is.data.frame(x))
        sapply(x, sd, na.rm = na.rm)
    else sqrt(var(as.vector(x), na.rm = na.rm))
}

とするだけで，どこがわるいのだろうか？？？

sd(データフレーム) は廃止

sd 関数は，もはや，行列やデータフレームを対象としない。

> class(d)
[1] "matrix"
> sd(d)
[1] 1 1
 警告メッセージ：
sd() is deprecated.
 Use apply(*, 2, sd) instead.
> sd(data.frame(d))
X1 X2
 1  1
 警告メッセージ：
sd() is deprecated.
 Use sapply(*, sd) instead.

なぜ，結果が character になるのだろうか

> x <- data.frame(x=1:5, y=c(3,2,5,4,3))
> sapply(x, factor)
     x   y 
[1,] "1" "3"
[2,] "2" "2"
[3,] "3" "5"
[4,] "4" "4"
[5,] "5" "3"

cxxfunction 最も単純な（？）平均値を求める場合

平均値を求めるくらいなら，何の差もないだろうと...

最初に，素直にベクトル要素に添え字でアクセスしてみたら，mean 関数の速度の半分ほど（要するに，遅いと言うこと）。

「これではならじ」と，ポインタを使って書いたら，

src <- '
    Rcpp::NumericVector x(X);
    int n = x.length();
    double *p;
    double sum = 0;
    p = REAL(x);
    for (int i = 0; i < n; i++) {
        sum += *p++;
    }
    return wrap(sum/n);
'
library(inline)
Mean <- cxxfunction(signature(X="numeric"), src, plugin="Rcpp")

倍速にはなった。

しかし，この場合も他と比べてなんという発見もない。

> x <- rnorm(100000000)

> system.time(a <- mean(x))
                  ユーザ                  システム 
0.38000000000000966338121 0.00099999999999944577667
                    経過 
0.38300000000003819877747
> system.time(b <- Mean(x))
               ユーザ               システム 
0.17300000000000181899 0.00000000000000000000
                 経過 
0.17199999999991177901

> dump(a)
a = -0.000169827373561765

> dump(b)
b = -0.000169827373561735

> a == b
[1] FALSE
> a-b
[1] -2.99781900692241976e-17

100000000 個のデータの平均値が，たかだか 0.4 秒なんだから。
シミュレーションで 1 日に約 20 万回やれるということです。

cxxfunction n! を全てリストアップする

e1071 に permutations という関数がある。
実装はいろいろあるけど，簡単そうな奥村先生に典拠するアルゴリズムを使って書いてみた。
ただし，n! なんて，すぐにメモリ制限を超えてしまうので，実際の利用はたかだか n <= 10 しかないので，どんなに速いアルゴリズムでも，生涯期待計算時間（ともいうようなもの）を見ると，ゼロに近いだろう。（いろいろなプログラムの実装法を経験するという意味はあるだろうとやってみているだけなのだが...それにしても，本当に利用すべきと言う場合は未だ見つからず）。

src <- '
    int N = as(n);
    if (N < 1 || N > 10) {
        return wrap(NA_REAL);
    }
    int i, k, t, c[N + 1], p[N];
    int size = 1;
    for (i = 2; i <= N; i++) {
        size *= i;
    }
    Rcpp::IntegerMatrix res(size, N);

    int count = 0;
    for (i = 0; i < N; i++) {
        p[i] = i + 1;
    }
    for (i = 1; i <= N; i++) {
        c[i] = i;
    }
    k = 1;
    while (k < N) {
        if (k & 1) {
            i = c[k];
        }
        else {
            i = 0;
        }
        t = p[k];
        p[k] = p[i];
        p[i] = t;
        for (i = 0; i < N; i++) {
            res(count, i) = p[i];
        }
        count++;
        k = 1;
        while (c[k] == 0) {
            c[k] = k;
            k++;
        }
        c[k]--;
    }
    return wrap(res);
'
library(inline)
system.time(genperm <- cxxfunction(signature(n="integer"), src, plugin="Rcpp"))
system.time(genperm(10))
library(e1071)
system.time(permutations(10))

実行例

> system.time(genperm <- cxxfunction(signature(n="integer"), src, plugin="Rcpp"))
   ユーザ   システム       経過 
     3.146      0.708      8.395 # コンパイルするのに必要な時間

> system.time(genperm(10)) # 10! のリストを計算する時間
   ユーザ   システム       経過 
     0.253      0.123      0.373

> library(e1071)
 要求されたパッケージ class をロード中です
> system.time(permutations(10)) # perumutations による 10! のリストを計算する時間
   ユーザ   システム       経過 
     1.537      0.589      2.109

ほとんど無駄な努力と言える。

ダメ出し：2011-11-06 ２次計画問題

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)

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λ

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

cxxfunction NA の扱いと，ポインターの利用

移動平均は 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

ダメ出し：hatena グループの閉鎖性

コメントするのにログインまたは登録しろと。。。
そんなのやだね。登録者数を増やしたいのか。

http://d.hatena.ne.jp/tsutatsutatsuta/20111024

行列の各要素ごとに異なる引数を与えて，同じ関数を適用

applyを使って，行列の行や列ごとに，異なる引数を与えて同じ関数を適用したいときの方法です．

まずこんな行列を作成

> (xx <- t(array(1:12, dim = c(4, 3))))
     [,1] [,2] [,3] [,4]
[1,]    1    2    3    4
[2,]    5    6    7    8
[3,]    9   10   11   12
> 

列ごとに異なる引数を与えて，各行に同じ関数を適用したいとき．

例えば，第1列なら1，第2列なら2，第3列なら3，第4列なら4を引きたいとき．

> t(apply(xx, 1, function(x){x - 1:4}))
     [,1] [,2] [,3] [,4]
[1,]    0    0    0    0
[2,]    4    4    4    4
[3,]    8    8    8    8
> 

行ごとに異なる引数を与えて，各列に同じ関数を適用したいとき．

例えば，第1行なら0，第2行なら4，第3行なら8を引きたいとき．

> apply(xx, 2, function(x){x - c(0, 4, 8)})
     [,1] [,2] [,3] [,4]
[1,]    1    2    3    4
[2,]    1    2    3    4
[3,]    1    2    3    4
> 

それだけです．
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

それだけです？

そんな大げさなことする必要はさらさらない。以下で十分

> t(t(xx)-1:4)
     [,1] [,2] [,3] [,4]
[1,]    0    0    0    0
[2,]    4    4    4    4
[3,]    8    8    8    8

> xx-c(0, 4, 8)
     [,1] [,2] [,3] [,4]
[1,]    1    2    3    4
[2,]    1    2    3    4
[3,]    1    2    3    4

> xx-0:2*4
     [,1] [,2] [,3] [,4]
[1,]    1    2    3    4
[2,]    1    2    3    4
[3,]    1    2    3    4