オレオレ平滑化をやってみた

この前自前にて密度トレイスを書いてみましたが、これを使えば平滑化も簡単に行けるということが分かったので、実際にやってみました。

ガウシアンカーネル関数はこの前とほとんど同じで、 $\int y f_{X,Y}(x,y)dy$ を計算するための関数を作りました。あとはそれを割ったりしているだけ。簡単にできました。

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

y.exp <- function(X,Y,h){
  y.exp.intern <- function(X,Y,h){
    return(paste("1/sqrt(2 * pi) * exp(-1/2 * ((x - ",X,")/h)^2) *",Y,sep=""))
  }
  eval(parse(text=paste("function(x){1/h * ",paste(mapply(function(X,Y){y.exp.intern(X,Y,h)},X,Y),collapse="+"),"}")))
}

x <- cars$speed
y <- cars$dist
h <- 1
plot(x,y,xlab="speed of cars",ylab="speed of dist")
title("Rの平滑化関数と自前の関数の比較")
lines(lowess(x,y),col="red")
s <- seq(min(x),max(x),length.out=1000)
g <- sapply(s,gauusian(x,h))
y <- sapply(s,y.exp(x,y,h))
lines(s,y/g,col="blue")
legend(5,100,c("Rのlowess","自前の平滑化曲線"),lwd=1,col=c("red","blue"))

f:id:syou6162:20161024230323p:plain

この場合だとえらく自前とRのlowessが違うように見えるのですが、バンド幅を変えてやれば似たような感じになります。バンド幅を2倍にしてみました。
f:id:syou6162:20161024230341p:plain

色々な基準を計算して、最適なバンド幅をディフォルトに指定しておけばよいわけですが、それはまた今度。Rのlowessとかだと2つのベクトルしか返してくれないので、信頼区間とかの計算があれになってくるわけですが、自前でやれるようにしておくと色々できるのであとで役に…立つかもしれない。

追記

上のは関数を文字列として出しまくって、evalしているんだけどもっと簡単にしてみた。たぶん、こっちのほうが早い気がする。

こっちは関数が要素のものをlistに、sapplyでそれぞれの関数ごとに値を出す→applyでsumを使って、各xの値ごとに予測値を出すということをやっている。

gauusian <- function(X,h){
  return(function(x){1/(length(X)*h) * 1/sqrt(2 * pi) * exp(-1/2 * ((x - X)/h)^2)})
}

y.exp <- function(X,Y,h){
  return(function(x){1/(length(X)*h) * 1/sqrt(2 * pi) * exp(-1/2 * ((x - X)/h)^2) * Y})
}

x <- cars$speed
y <- cars$dist
h <- 0.5
plot(x,y,xlab="speed of cars",ylab="speed of dist")
title("Rの平滑化関数と自前の関数の比較")
lines(lowess(x,y),col="red")
s <- seq(min(x),max(x),length.out=1000)
g <- apply(sapply(s,gauusian(x,h)),2,sum)
y <- apply(sapply(s,y.exp(x,y,h)),2,sum)
lines(s,y/g,col="blue")
legend(5,100,c("Rのlowess","自前の平滑化曲線"),lwd=1,col=c("red","blue"))

密度トレイスのやつも一応書いておく。

gauusian <- function(X,h){
  return(function(x){1/(length(X)*h) * 1/sqrt(2 * pi) * exp(-1/2 * ((x - X)/h)^2)})
}
s <- seq(min(x)-sd(x),max(x)+sd(x),length.out=1000)
plot(s,apply(sapply(s,gauusian(x,2.15)),2,sum),type="l")