昨日の続きのところ。Cでナダラヤワトン推定量*1計算させて、Rのoptimを使って計算したら早くなったっぽい気がした。でも、一変数はRだけでまともなコードがなかった(文字列evalとかのやつしか残っていなかった)ので速度比較ができなかった。
これじゃ最適化関数もCで書かないといけないのか、そこまでしなくていいいのかの判断材料がないので、Rで書いて比較することにした。Cのプログラムは昨日のところのまま。
dyn.load("/tmp/hoge.so") cv <- function(y,x1,h) { .C("cv", as.double(y), as.double(x1), as.integer(length(y)), as.double(h), result = 0)$result } m <- function(y,x,x1,h) { .C("m", as.double(y), as.double(x), as.double(x1), as.integer(length(y)), as.double(h), result = 0)$result } epanechnikov <- function(X1,h1){ X1 <- X1 return(function(x1){ u <- (x1 - X1)/h1 if(1 > abs(u)){ return(0.75 * (1-u^2)) }else{ return(0) } }) } y.exp <- function(X1,Y,h1){ X1 <- X1 Y <- Y return(function(x1){ u <- (x1 - X1)/h1 if(1 > abs(u)){ return(0.75 * (1-u^2) * Y) }else{ return(0) } }) } n <- 1000 s <- seq(-30,30,length.out=n) y <- (function(x){sin(x)})(s) + rnorm(n,0,0.1) plot(s,y) (h <- optimize(function(h){cv(y,s,h)},c(0.1,0.4))$minimum) lines(s,mapply(function(x){m(y,s,x,h)},s)) lines(lowess(s,y),col="red") mean.h <- function(n){ X1 <- seq(-3,3,length.out=n) Y <- (function(x){sin(x)})(X1) + rnorm(n,0,0.1) min <- 0.001 max <- 1 t_c <- system.time(h_c <- optimize(function(h){cv(Y,X1,h)},c(min,max))$minimum)[3] t_r <- system.time(h_r <- optimize(function(h1){ sum(mapply(function(X1i,Yi){ (Yi - (Reduce("+",Map(function(f){f(X1i)},mapply(function(X1,Y){y.exp(X1,Y,h1)},X1,Y))) - 0.75*Yi) / (Reduce("+",Map(function(f){f(X1i)},mapply(function(X1){epanechnikov(X1,h1)},X1))) - 0.75) )^2 },X1,Y)) },c(min,max))$minimum)[3] cat(h_c,"\t",h_r,fill=TRUE) time <- c(t_c,t_r) names(time) <- NULL return(time) }
計算時間の短縮
プログラム間違ってるんじゃないかと思ったので、Cでのクロスバリデーションを最小にするhとRでのhとを出力させている。5回の計算時間の平均を求めるとこんな感じに。> apply(mapply(mean.h,rep(100,5)),1,mean) 0.1893555 0.1893555 0.2037649 0.2037649 0.2211508 0.2211508 0.1818008 0.1818008 0.2865864 0.2865864 [1] 0.0040 5.5206 > apply(mapply(mean.h,rep(200,5)),1,mean) 0.1615935 0.1615935 0.1903541 0.1903541 0.3099866 0.3099866 0.3015350 0.3015350 0.2411889 0.2411889 [1] 0.0134 22.2692 > apply(mapply(mean.h,rep(400,5)),1,mean) 0.2658077 0.2658077 0.2255306 0.2255306 0.2060307 0.2060307 0.2256248 0.2256248 0.2548975 0.2548975 [1] 0.0572 106.5370 > apply(mapply(mean.h,rep(1000,5)),1,mean) 0.2054734 0.2054734 0.1681855 0.1681855 0.1848514 0.1848514 0.1892796 0.1892796 0.1982217 0.1982217 [1] 0.3398 769.6374
| サンプル数 | n=100 | n=200 | n=400 | n=1000 |
|---|---|---|---|---|
| Cを使った場合の計算時間 | 0.0040 | 0.0134 | 0.0572 | 0.3398 |
| Rを使った場合の計算時間 | 5.5206 | 22.2692 | 106.5370 | 769.6374 |
| 計算時間比(R/C) | 1380倍 | 1661倍 | 1862倍 | 2264倍 |
もちろん同じデータを使っているので、最適解も同じになるはずで上のように同じになっている。そんでもって驚くべきは計算時間がアホみたいに早くなっているということだ。1000倍以上早くなっているし、データ数が多くなるほど、比が大きくなっている。
なんかうまく行きすぎてて、たぶんどっかでミスってると思うんだけどどうなのか。最適化をするための関数を自分で用意するのは非常に手間がかかる。実装するとなると、最適化手法は初期点の問題や収束性の問題があるから一つの手法だけではなく、いくつか作って試すべきだ(最初は初期点に依存しないような手法であたりを付けて、それを元に鋭い刃物を使う、というような)。そうなると、勉強しないといけないことが増えまくるし、自分が作ったものはバグが入る可能性が大きい(笑…ちゃいけないけど)。
最適化してあげたい関数の中身だけCで書いてあげて(これは昨日のコードを見ると分かるけど、ものすごーく簡単)、信頼のおけるRの最適化関数にあとは投げられて、しかもRだけの時よりものすごく早い。こんなおいしい話があるときには注意しなさいよ、と小学校の頃から言われているんだけど、どうなのか。コンパイルするとそんなに早くなるものなのか?とりあえず色んな人に聞いてみよう。
*1:leave-one outしてるので、正確には違うのだが
