カーネル平滑化のメモ

今回はカーネル平滑化による曲線フィッティングについてである。カーネル平滑化というのは、簡単に言うと、x軸の各値で、周囲の点にガウス分布の形をした重みをかけて、y軸の値を推定するノンパラメトリックな平滑化手法である。

重みというのは、例えばx=3.2の点を見ており、周囲にx=1,3,4のデータがあれば、データに対する重みは x=3のデータ > x=4のデータ > x=1のデータの順にするということである。ガウス分布の裾を大きくすると、x軸のどの値でもデータ全体を見るようになるので、次第にx軸に平行な平滑化の直線が得られることが予想される。KernSmoothというライブラリのksmoothでカーネル平滑化は実行できるので、example(ksmooth)に手を加えて、それを確かめてみよう。

例に用いられているデータcarsは、車の速度と車が減速してから停止するまでに移動した距離の関係を表したものである。

library(KernSmooth)
require(graphics)

colset <- 2:6
bandset <- c(2,5,10,20,100)

# for文を書かないで1要素ずつ処理
myfunc <- function(bc,xx,yy,type="normal"){
 pp <- ksmooth(xx,yy,type,bandwidth=bc[1])
 lines(pp,,col=bc[2],lwd=2)
 return(pp)
}

plot(cars$speed,cars$dist);
targetY <- 50
abline(h=targetY)
bc <- cbind(bandset,colset)
coords <- apply(bc,1,myfunc,cars$speed,cars$dist)

legend(x=min(cars$speed),y=max(cars$dist),
 legend=paste("bd",bandset,sep=""),pch=19,col=colset)

以上の内容がフィッティングである。次に、このフィッティングカーブを持ちいて、距離から速度を推定することを行う。

ksmooth関数はx軸とy軸の値をセットにして返すので、lines関数で曲線に近い直線を描くことができる。その精度は、オプションで、点数を大きい値に指定することによって実現できる。

y軸の値に最も近いx軸のインデックスを返す関数を以下に示す。実行するときには、上のコードの下に付け加えればよい。

# y軸の値がvalのときのx軸の値を返す
# 初めて超えた点を近似値として返す
retEstIntersection <- function(vec,val){
 tarIdx <- which(vec$y>val)
 if( length(tarIdx) < 1 ){
  return(0)
 }
 idx <- min(which(vec$y>val))
 return(idx)
 # return(which.min(abs(vec$y-val))) # 一番近い点を求める場合
}

# targetYに最も近いデータの行番号が返ってくる
# これをksmoothの返り値のリストに放り込めば，交点のx軸の近似値が返せる
ret <- sapply(coords,retEstIntersection,targetY,simplify=T)
# 上の行はこんな書き方でも代用できる
# ret <- lapply(coords,retEstIntersection,targetY,simplify=T)
# result <- unlist(ret)

以上のコードを実行すると、下記のものが得られる。

1)y軸の値がtargetYを初めて超えたx軸のインデックスを返した結果

> ret
[1] 64 67 70 80  0

ここで0は平滑化で得られた曲線とy=targetYと交わらなかったことを意味する。

2)プロット

フィッティングの精度を見てみると、bandwidthが小さいとデータに引っ張られやすく、大きいとx軸に平行な直線に近いものになり、精度が悪くなる。今回の例では、bandwidth=5もしくは10が一番きれいにフィッティングできているように見える。