クロスバリデーション最小化のコード関係をCで書いてみた

Cで書いてどうなるかを見てみるテスト。やってみて気がついたが、サンプル数が少ないとき、epanechnikovだと分母が0の領域が出てきて微妙な感じになってしまうのでどうにかしないといけないが、面倒なので放置している。。。サンプル少ない時はガウシアンカーネルで逃げることにするか(ぉ。

あとCで書いたらどれくらい早くなるのかを試そうと思った、ということで書いているのだが、一変数の時の最適化はまともなコードじゃないので一概に比較できないという感じである。

carのデータをナダラヤワトソン推定量で回帰させるのと、sinからのデータのjitteredしたデータを回帰、という二種類をやっている。

#include <R.h> 
#include <Rinternals.h> 

double epanechnikov(double u)
{
  if(-1.0 < u && u < 1.0){
    return ((3.0 / 4.0) * (1.0 - u*u));
  }else{
    return 0;
  }
}

void m(double *y,double *x1,double *x,int *n,double *h,double *result)
{
  double numerator = 0,denominator = 0;
  for(int i=0;i<*n;i++){
    numerator += epanechnikov(((*x - x1[i]) / *h)) * y[i];
  }
  for(int i=0;i<*n;i++){
    denominator += epanechnikov((*x - x1[i]) / *h);
  }
  *result = numerator / denominator;
}

double g_i(double *y,double *x1,double x,int n,double h,int i)
{
  double sum = 0;
  for(int j=0;j<n;j++){
    if(j!=i) sum += epanechnikov(((x - x1[j]) / h)) * y[j];
  }
  return sum;
}

double f_i(double *x1,double x,int n,double h,int i)
{
  double sum = 0;
  for(int j=0;j<n;j++){
    if(j!=i) sum += epanechnikov((x - x1[j]) / h);
  }
  return sum;
}

double m_i(double *y,double *x1,double x,int n,double h,int i)
{
  return g_i(y,x1,x,n,h,i) / f_i(x1,x,n,h,i);
}

void cv(double *y,double *x1,int *n,double *h,double *result)
{
  int i,j;
  for(i =0;i<*n;i++){
    *result = *result + pow(y[i] - m_i(y,x1,x1[i],*n,*h,i),2);
  }
}

次がR側のコード。一変数最適化はoptimじゃなくてoptimizeを使えと怒られた。

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 
}

attach(cars)
s <- seq(5,25,length.out=1000)
plot(speed,dist)
h <- optimize(function(h){cv(dist,speed,h)},c(1,10))$minimum
lines(s,mapply(function(x){m(dist,speed,x,h)},s))
lines(lowess(speed,dist),col="red")

n <- 200
s <- seq(-3,3,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){
  s <- seq(-3,3,length.out=n)
  y <- (function(x){sin(x)})(s) + rnorm(n,0,0.1)
  (h <- optimize(function(h){cv(y,s,h)},c(0.1,0.4))$minimum)
  return(h)
}

carsだとlowessとそんなに変わらないけど、sinだとlowessいまいちだなという感じか。overfittingしているとも言えなくもないが、基準次第か。


かつで、ならばナダラヤワトソン推定量は真の関数に確率収束する、などと言っているので、サンプル増やした時のhの挙動を見てみることにした。一回だとぶれてよく分からないので、それぞれのサンプル数において20回の平均を算出した。遅いけど、hは0に向かっているようだと言えそうである。

> mean(mapply(mean.h,rep(100,20)))
[1] 0.2979589
> mean(mapply(mean.h,rep(150,20)))
[1] 0.2761686
> mean(mapply(mean.h,rep(200,20)))
[1] 0.2647541
> mean(mapply(mean.h,rep(300,20)))
[1] 0.2255239
> mean(mapply(mean.h,rep(1000,20)))
[1] 0.1883825
> mean(mapply(mean.h,rep(5000,20)))
[1] 0.1232544

20も回すと面倒なので、2回だけど万単位で回すとこんな感じ。ふむ。

> mean(mapply(mean.h,rep(10000,2)))
[1] 0.1015914
> mean(mapply(mean.h,rep(20000,2)))
[1] 0.08564092
> mean(mapply(mean.h,rep(40000,2)))
[1] 0.06247027