Confidence interval for kernel density estimator with optimal bandwidth under the cross-validation criterion

ここの続き。カーネル密度推定の場合においての最適なバンド幅をクロスバリデーションで求める、というのはまだやっていなかったのでやりました。

カーネル密度推定の場合のクロスバリデーションは、積分の計算が入ってくるので、数値積分をやったりしました。Rにはintegrateという数値積分してくれる関数があるんですが、C側で何回も使いたかったので、自前で書いてみました。カーネル密度推定のクロスバリデーション関数は
\begin{align}\mbox{CV}(h) &= \frac{1}{n^2h^2} \sum_{i=1}^{n}\sum_{j=1}^{n} \int K(\frac{x-x_i}{h}) K(\frac{x-x_j}{h}) dx - \frac{2}{n(n-1)h} \sum_{i=1}^{n}\sum_{j=1,j \neq i}^{n} K(\frac{x_i-x_j}{h}) \\ &= \frac{1}{n^2h} \sum_{i=1}^{n}\sum_{j=1}^{n} K \bigcirc K(\frac{x_i-x_j}{h}) - \frac{2}{n(n-1)h} \sum_{i=1}^{n}\sum_{j=1,j \neq i}^{n} K(\frac{x_i-x_j}{h})\end{align}
という感じになってます。

とりあえずクロスバリデーションに必要になるCのコード。こいつをoptimize関数に投げます。

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

double epanechnikov(double X, double x,double h)
{
  double u = (x - X) / h;
  if(-1.0 < u && u < 1.0){
    return 3.0 / 4.0 * (1.0 - u * u);
  }else{
    return 0.0;
  }
}

double convolution(double Xi,double Xj,double x,double h)
{
  // 畳み込みの中身の関数
  double ui = (x - Xi) / h;
  double uj = (x - Xj) / h;
  if(-1.0 < ui && ui < 1.0){
    if (-1.0 < uj && uj < 1.0){
      return 3.0 / 4.0 * (1.0 - ui * ui) * 3.0 / 4.0 * (1.0 - uj * uj);
    }
  }else{
    return 0.0;
  }
}

double integrate(double Xi,double Xj,double h,int n)
{
  // 数値積分をする関数
  int i;
  double a,b;

  double ak[n];
  if(Xi - h < Xj - h){
    a = Xi - h;
  }else{
    a = Xj - h;
  }

  if(Xi + h > Xj + h){
    b = Xi + h;
  }else{
    b = Xj + h;
  }

  double diff = (b - a) / n;
  double sum = a;

  double conv_i = 0.0;
  double conv_i_minus1 = 0.0;

  for (i = 0; i < n; i++){
    ak[i] = sum;
    sum += diff;
  }
  ak[n-1] = sum;
 
  sum = 0.0;
  for (i = 1; i < n; i++){
    conv_i_minus1 =   convolution(Xi,Xj,ak[i-1],h);
    conv_i        =   convolution(Xi,Xj,ak[i],h);
    sum += (ak[i] - ak[i-1]) * (conv_i_minus1 + conv_i) / 2.0;
  }

  return sum;
}

void cv(double *X,double *h,int *precision,int *n,double *result)
{
  // precisionは数値積分の区間の幅。狭いほうが精度がいい
  int i,j;
  double sum_part_one = 0.0;
  double sum_part_two = 0.0;
  for (i = 0; i < *n; i++){
    for (j = 0; j < *n; j++){
      sum_part_one += integrate(X[i],X[j],*h,*precision);
    }
  }
  // 変数変換する前のものを考えているので、n^2h^2で割る、というものでよい
  sum_part_one = sum_part_one / (*n * *n * *h * *h);

  for (i = 0; i < *n; i++){
    for (j = 0; j < *n; j++){
      if (j != i){
	sum_part_two += epanechnikov(X[i],X[j],*h);
      }
    }
  }
  sum_part_two = 2.0 / (*n * (*n - 1) * *h) * sum_part_two;

  *result =  sum_part_one - sum_part_two;
}

R側は前回と大体同じだけど、最適化の部分が入ってきてます。一変数の最適化の場合にはoptimじゃなくてoptimizeがいいんだってさ。あと、クロスバリデーション関数自体もplotしてみたり。

epanechnikov <- function(X1,h1){
  X1 <- X1
  h1 <- h1
  return(function(x1){
    u <- (x1 - X1)/h1
    if(1 > abs(u)){
      return(0.75 * (1-u^2))
    }else{
      return(0)
    }
  })
}

density.estimator <- function(X1,h,length.out=30){
  s <- seq(from=min(X1),to=max(X1),length.out=length.out)
  1 / (length(X1) * h) * apply(sapply(s,
               function(x){
                 mapply(
                        function(f){
                          f(x)
                        },mapply(
                                 function(X1){
                                   epanechnikov(X1,h)
                                 },X1)
                        )
               }
               ),2,sum)
}

confidence.interval <- function(X1,h,length.out=30){
  s <- seq(from=min(X1),to=max(X1),length.out=length.out)
  fhatx <- density.estimator(X1,h=h,length.out=length.out)
  upper <- fhatx + 1.96 * (length(X1)*h)^(-1/2) * (fhatx * 3/5)^(1/2)
  lower <- fhatx - 1.96 * (length(X1)*h)^(-1/2) * (fhatx * 3/5)^(1/2)
  return(list(fhatx=fhatx,upper=upper,lower=lower))
}

dyn.load("/tmp/hoge.so")
cv <- function(x,h,precision=1000){
  .C("cv",
     as.double(x),
     as.double(h),
     as.integer(precision),
     as.integer(length(x)),
     result = 0)$result 
}


X <- rnorm(200)
s <- seq(0,3,length.out=100)
plot(s,sapply(s,function(h)cv(X,h,100)),type="l")

n <- 100
s <- seq(from=min(X),to=max(X),length.out=n)
h <- optimize(function(h){cv(X,h,100)},c(0.5,1.5))$minimum
c <- confidence.interval(X,h=h,length.out=n)
d <- c$fhatx
u <- c$upper
l <- c$lower
plot(s,d,ylim=c(min(d)-0.02,max(d)+0.05),type="l")
lines(s,u,col="red")
lines(s,l,col="red")

クロスバリデーション関数の外形はこんな感じ。これでとりあえず当たりを付けておいて、最適化の時の材料にすることにします。ナダラヤワトソンのときも同じようなことをやっていたり
Quartz 2 [*]
で、クロスバリデーション関数を最小にするバンド幅を使った時の信頼区間はこんな感じ。まあ、いいかな。
Quartz 2 [*]-11
最適化の処理をするときに、目的関数をCレベルで書いて高速化するっていうのはもう大分慣れてきたな。よいことだ。