カーネル密度推定メモ

# [-h,h]の範囲の一様分布をカーネル関数とする
my.rect <- function(x,h){
	ret <- rep(0,length(x))
	ret[which(abs(x)<=1)] <- 1
	ret
}
hs <- c(0.1,0.25,0.5,1)
y <- seq(from=-10,to=10,length=1000)
x.k <- matrix(0,length(y),length(hs))
for(i in 1:length(hs)){
	h <- hs[i]
	for(j in 1:length(x)){
		x.k[,i] <- x.k[,i] + 1/(N*h)*my.rect((y-x[j])/h,h)
	}
}
matplot(y,x.k,type="l")
legend(-9,0.5,hs,lty=1:length(hs),col=1:length(hs))

x <- cbind(rnorm(100,0,1),rnorm(100,0,2))
x <- rbind(x,cbind(rnorm(50,3,1),rnorm(50,5,1)))
plot(x)
@
<<fig=TRUE>>=
h <- 0.5
X <- seq(from=-5,to=5,length=100)
Y <- matrix(0,length(X),length(X))
for(k in 1:length(x)){
	# ２次元の軸ごとにカーネル分布を算出しそれを掛けあわせてある
	tmp <- outer(dnorm((x[k]-X)/h),dnorm((x[k]-X)/h),"*")
	Y <- Y + tmp
}
image(Y)
contour(Y,add=TRUE)