多変数における最小二乗クロスバリデーション評価関数を最小にするhを決定する

以前カーネル単回帰のところで最小二乗クロスバリデーション評価関数
$CV(h_x) = \sum^n_{i=1}(Y_i - (\frac{\sum^n_{j=1,j \neq i}K_X(\frac{X_i-X_j}{h_x})Y_j}{\sum^n_{j=1,j \neq i} K_X(\frac{X_i-X_j}{h_x})}))^2$
を最小にするhを決めるというのをやりました。

「カーネル重回帰でも最小二乗クロスバリデーション評価関数を最小にするのがある意味(hogehogeに確率収束する、など)においてよい」というのはないんですが*1、まあ何かしらの多変数関数を最小にしないといけないので、練習がてら多変数における最小二乗クロスバリデーション評価関数の最小化問題に取り組むことにしました。

乗法的核関数
それぞれのカーネルはガウシアンカーネル
とする
- 要するにp個の変数でやるということ

という風にすると、最小二乗クロスバリデーション評価関数は以下のように書けます。
$\begin{eqnarray}CV(H) &=& \sum^n_{i=1}(Y_i - (\frac{\sum^n_{j=1,j \neq i}K(H^{-1}(X_i-X_j))Y_j}{\sum^n_{j=1,j \neq i} K(H^{-1}(X_i-X_j))}))^2 \\ &=& \sum^n_{i=1}[Y_i - (\frac{\sum^n_{j=1,j \neq i}\prod^p_{k=1}{\exp(-\frac{1}{2}\{\frac{(X_{ik}-X_{jk})}{h_k}\}^2)}Y_j}{\sum^n_{j=1,j \neq i}\prod^p_{k=1}{\exp(-\frac{1}{2}\{\frac{(X_{ik}-X_{jk})}{h_k}\}^2)}})]^2 \end{eqnarray}$

前のコードでは、

最小二乗クロスバリデーション評価関数を文字列として生成
それをD関数で数式微分
evalなどをして、最急降下法に持っていく

という流れでやっていました。しかし、

関数が巨大すぎて、文字列の連結に結構時間がかかる
- 一変数でこれなので、多変数でそのままつっぱしると危険
最急降下法以外の最適化関数もちょっと使っておきたい

ということで、基本的に「文字列化は使わないで、Rに付属の最適化関数を使う」という方針で行くことにしました。そういうわけで以下がそのコード。sum以下が(長いですが)全部関数です。mapply(に相当するもの)が3つくらいネストしたりしていますが、怖がってはいけません。

l <- 10
X1 <- seq(0,1,length=l)
X2 <- seq(0,1,length=l)

e <- expand.grid(X1,X2)
Y <- apply(e,1,function(x){
  x1 <- x[1]
  x2 <- x[2]
  mkr_with_noise(x1,x2)
})
X1 <- e[,1]
X2 <- e[,2]

optim(par=c(1,0.1),function(h){
  h1 <- h[1]
  h2 <- h[2]
  sum(mapply(function(X1i,X2i,Yi){
    (Yi -
     (Reduce("+",Map(function(f){outer(X1i,X2i,f)},mapply(function(X1,X2,Y){y.exp(X1,X2,Y,h1,h2)},X1,X2,Y))) - 1/sqrt(2*pi)) /
     (Reduce("+",Map(function(f){outer(X1i,X2i,f)},mapply(function(X1,X2){gauusian(X1,X2,h1,h2)},X1,X2))) - 1/sqrt(2*pi))
     )^2
  },X1,X2,Y))
})

これを動かすと、こんな感じで最適なhのベクトルが返ってきます。

> optim(par=c(0.1,0.1),function(h){
+   h1 <- h[1]
+   h2 <- h[2]
+   sum(mapply(function(X1i,X2i,Yi){
+     (Yi -
+      (Reduce("+",Map(function(f){outer(X1i,X2i,f)},mapply(function(X1,X2,Y){y.exp(X1,X2,Y,h1,h2)},X1,X2,Y))) - 1/sqrt(2*pi)) /
+      (Reduce("+",Map(function(f){outer(X1i,X2i,f)},mapply(function(X1,X2){gauusian(X1,X2,h1,h2)},X1,X2))) - 1/sqrt(2*pi))
+      )^2
+   },X1,X2,Y))
+ })
$par
[1] 0.06720100 0.08223122

$value
[1] 4.481639

$counts
function gradient 
      57       NA 

$convergence
[1] 0

$message
NULL

optim関数はディフォルトでNelder-Mead法というの使っているようです。Nelder-Mead法というのは、

関数値だけを使う
頑健(例えば初期値の選択に敏感でない)であるが、相対的に遅い
微分できない関数に対してもそれなりに使える

などの特徴があるそうです。

で、このhを使って多変数ナダラヤワトソン推定量をplotしてみました。まあ、悪くはないかなーという感じかな。幅が荒いのは、細かくすると最適化に時間が結構かかるからやってないだけですwww。

以下そのコード。

l <- 10
X1 <- seq(0,1,length=l)
X2 <- seq(0,1,length=l)
Y <- outer(X1,X2,mkr_with_noise)

e <- expand.grid(X1,X2)
Y <- apply(e,1,function(x){
  x1 <- x[1]
  x2 <- x[2]
  mkr_with_noise(x1,x2)
})
X1 <- e[,1]
X2 <- e[,2]

h1 <- 0.06720100 
h2 <- 0.08223122

s1 <- seq(0,1,length=l)
s2 <- seq(0,1,length=l)

persp(s1,s2,
      Reduce("+",Map(function(f){outer(s1,s2,f)},mapply(function(X1,X2,Y){y.exp(X1,X2,Y,h1,h2)},X1,X2,Y))) /
      Reduce("+",Map(function(f){outer(s1,s2,f)},mapply(function(X1,X2){gauusian(X1,X2,h1,h2)},X1,X2))),
      main=paste("多変量ナダラヤワトソン推定量(h1=",h1,",h2=",h2,")\n",
        "(乗法的カーネル関数を仮定し、h1はCV関数の最小化によって決定)"),
      theta=320,phi=20,col=rainbow(50),ticktype="detailed",
      xlab="x1",ylab="x2",zlab="density")

*1:あるかもだけど、僕は今のところ知らない