カーネル関数を使って密度トレイスを書いてみた→最後のほうに微妙に間違いがあるので注意!!

カーネル関数の一つであるエパネクニコフというやつをとりあえずやってみた。Rでいうところのdensityを自前で書いた、という感じです。関数は6行なんだけど、すっごい重い。

epanechnikov <- function(r){
  epanechnikov.intern <- function(r){
    return(paste("3/4 * (1 - (x - ",r,")^2) * ifelse(((x < 1.0 + ",r,") && (x > -1.0 + ",r,")),1,0)",sep=""))
  }
  eval(parse(text=paste("function(x){",paste(sapply(r,epanechnikov.intern),collapse="+"),"}")))
}

r <- rnorm(30)
s <- seq(-3,3,length.out=1000)
plot(s,sapply(s,epanechnikov(r)),type="l")

y軸の値が密度のやつになってないや。どうにかしよう。

追記

カーネル関数も密度関数であり、全区間で積分すれば1になる。だから、カーネル関数の和の積分はデータ数と一致する、ということに気がつけば簡単だった。

plot(density(r))
r <- rnorm(30)
s <- seq(-3,3,length.out=1000)
e <- sapply(s,epanechnikov(r))
plot(s,e/length(r),type="l")

ついでにガウシアンカーネル関数も書いてみた。エパネクニコフってのはifelseがデータ数分あって重かったけど、条件分岐がいらないからかこっちだと割りと軽い感じだった。

gauusian <- function(r){
  gauusian.intern <- function(r){
    return(paste("1/sqrt(2 * pi) * exp(-1/2 * (x - ",r,")^2)",sep=""))
  }
  eval(parse(text=paste("function(x){",paste(sapply(r,gauusian.intern),collapse="+"),"}")))
}

r <- rnorm(100,0,10)
s <- seq(min(r),max(r),length.out=1000)
e <- sapply(s,gauusian(r))
plot(s,e/length(r),type="l")

追記の追記

嘘、ごめん、これ嘘だわwww。いや、正確には嘘とまでいかないんだけど、違うという感じ。supportの値を固定してしまっていたようだ。勝手にh=1と勘違いしていた。

ノンパラだとこのsupportをいくつにするかが結構大きな問題っぽいから勝手に1に固定にしてたら発展できない。Rのdensityもそうだったことを思い出すんだ!!というわけで速攻訂正版。supportを0.1とかにするとはりねずみになったwww。

gauusian <- function(r,h){
  gauusian.intern <- function(r,h){
    return(paste("1/sqrt(2 * pi) * exp(-1/2 * ((x - ",r,")/h)^2)",sep=""))
  }
  eval(parse(text=paste("function(x){1/h * ",paste(sapply(r,function(s){gauusian.intern(s,h)}),collapse="+"),"}")))
}
r <- rnorm(1000,0,10)
s <- seq(min(r),max(r),length.out=1000)
e <- sapply(s,gauusian(r,0.1))
plot(s,e/length(r),type="l")