カーネル平滑化のメモ

今回はカーネル平滑化による曲線フィッティングについてである。カーネル平滑化というのは、簡単に言うと、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が一番きれいにフィッティングできているように見える。

applyに渡す関数を自前で書いてみる

Rの高速化のポイントにあげられるものとしてforを使わないということがある。
ベクトル単位で一度に処理せよということである。
今回はベクトル単位で処理する際によく用いられるapplyの挙動について紹介してみる。

まずは肩慣らしに、普通のapplyを使ってみよう。

> mat <- matrix(1:12,nrow=4)
> mat
     [,1] [,2] [,3]
[1,]    1    5    9
[2,]    2    6   10
[3,]    3    7   11
[4,]    4    8   12
> mat[2,3]<- NA
> mat
     [,1] [,2] [,3]
[1,]    1    5    9
[2,]    2    6   NA
[3,]    3    7   11
[4,]    4    8   12
> apply(mat,1,mean)
[1]  5 NA  7  8
> apply(mat,1,mean,na.rm=T)
[1] 5 4 7 8

関数にオプション(na.rm=T)を与えるときには、関数の後ろにオプションを書けばよい。
次に自前の関数を用意してapplyがどのような挙動になるか見てみる。

今回は、あるベクトルと行列に対し、完全に内容が一致する行のインデックスを
取り出す操作を行うことにしよう。

mat <-matrix(1:12,nrow=4)
mat[2,3] <- NA
pattern <- c(3,7,11) # 調べようとする内容
myfunc <- function(m,pat){
 print(c(is.matrix(m),is.vector(m)))
 print(m)
 return( which(all(m==pat)) )
}
apply(mat,1,myfunc,pattern)

自作の関数(myfunc)では、最初の引数mにapplyで一番最初の引数、
2番目の引数patには、FUNの後にある引数が対応している。
上記のスクリプトを実行すると、以下の結果が得られる。

> apply(mat,1,myfunc,pattern)
[1] FALSE  TRUE # print(c(is.matrix(m),is.vector(m)))の内容
[1] 1 5 9 # print(m)の内容
[1] FALSE  TRUE
[1]  2  6 NA
[1] FALSE  TRUE
[1]  3  7 11
[1] FALSE  TRUE
[1]  4  8 12
[[1]] # applyの返り値、integer(0)は長さ0の整数型ベクトル
integer(0)

[[2]]
integer(0)

[[3]] # マッチしたので、数値が返ってくる。しかし、3ではなく1になる
[1] 1

[[4]]
integer(0)

1が返ってきたという上記の内容から、実は1行ずつ判定しているということがわかる。
しかも関数の中で、最初の引数であった行列はベクトルにdropしているのである。
だから、3番目が一致しているということを判別したければ、さらに上記のコードに

拡張を行う必要がある。

res <- apply(mat,1,myfunc,pattern) # 上のコードの最終行
ret <- lapply(res,function(x) length(x)>0) # 対応する行は長さが1以上
idx <- which(unlist(ret))

とすると、idxには3が入って、めでたく必要な情報が得られたことになる。
今回はマトリックスという列の属性がすべて共通なものであったため、
容易に比較ができたが、データフレーム（列ごとに属性が違う）の場合は、
比較対象の属性の中で一番表現範囲の広い型にしてしまってから比較を
行わなければならない。
大雑把ではあるが、文字列型>実数型>整数型>論理型の順で表現できる
範囲が広くなると考えてよい。 入力がどうなるか判断がつかない場合は、
とりあえず文字列型に変換すればよい。 たとえばmatや調べたいパターンが
データフレームであった場合、以下のような修正が考えられるであろう。

myfunc <- function(m,pat){
 mm <- as.character(unlist(m))
 patpat <- as.character(unlist(pat))
 return( which(all(mm==patpat)) )
}

追記：上ではapplyの理解のためにmyfuncの返り値にwhichを付けているが、
実はつけない方が楽にかける。以下のように書けば、データ長で、
行番号の有無を判定する作業が無くて済む。

myfunc <- function(m,pat){
 return( all(m==pat) )
}
ret <- apply(mat,1,myfunc,pattern)
result <- which(ret)