以前多変数ナダラヤワトソン推定量をpltoしてみるとことプログラミムを書いたりしました。
この時は乗法的カーネル関数を仮定していたわけだけど、それを仮定しないで平滑化パラメータ行列を変化させると、多変数ナダラヤワトソン推定量をplotしたときにどういう変化があるかやってみな、と言われていたのでやってみました。
で、書いたのが下のコード。多次元でもいけるように書いたはいいんだけど、ものすごーく遅い。4枚の画像を吐くのに30分たっても終わらない。outerを使わないでforを回しているのがボトルネックとなっている。が、下のコードをouterに書き直しただけでは、動かない。関数がベクトルを投げても動く形になっていないからのようだ。
mkr_with_noise<- function(x1,x2){ return(sin(2*pi*x1)+x2+rnorm(1,0,1/4)) } epanechnikov.radial <- function(X,H){ X <- X H <- H return(function(x){ u <- as.vector(solve(H) %*% (x - X)) utu <- sqrt(sum(u^2)) ((1-utu) * ifelse(utu <= 1,1,0)) }) } y.exp <- function(X,Y,H){ X <- X Y <- Y H <- H return(function(x){ u <- as.vector(solve(H) %*% (x - X)) utu <- sqrt(sum(u^2)) ((1-utu) * ifelse(utu <= 1,1,0)) * Y }) } plot.3d <- function(s1,s2,n,H){ persp(s1,s2, Reduce("+",Map(function(f){ z <- matrix(1,n,n) for(x in 1:n){ for(y in 1:n){ z[x,y] <- f(c(s1[x],s2[y])) } } return(z) },mapply(function(X,Y){y.exp(X,Y,H)},X,Y))) / Reduce("+",Map(function(f){ z <- matrix(1,n,n) for(x in 1:n){ for(y in 1:n){ z[x,y] <- f(c(s1[x],s2[y])) } } return(z) },mapply(function(X){epanechnikov.radial(X,H)},X))), theta=320,phi=20,col=rainbow(50),ticktype="detailed", xlab="x1",ylab="x2",zlab="density") } n <- 30 X1 <- seq(0,1,length=n) X2 <- X1 X <- expand.grid(X1,X2) Y <- apply(X,1,function(X){ x1 <- X[1] x2 <- X[2] mkr_with_noise(x1,x2) }) # convert matrix to list X <- unlist(apply(X,1,function(x){list(x)}),recursive=FALSE) s1 <- seq(0,1,length=n) s2 <- seq(0,1,length=n) setwd("~/tex") H <- matrix(c(0.2,0,0,0.2),2,2,byrow=TRUE) png("nadaraya.multi1.png") plot.3d(s1,s2,n,H) dev.off() H <- matrix(c(0.1,0,0,0.2),2,2,byrow=TRUE) png("nadaraya.multi2.png") plot.3d(s1,s2,n,H) dev.off() H <- matrix(c(0.1,0.2,0.2,0.1),2,2,byrow=TRUE) png("nadaraya.multi3.png") plot.3d(s1,s2,n,H) dev.off() H <- matrix(c(0.1,0.15,0.1,0.2),2,2,byrow=TRUE) png("nadaraya.multi4.png") plot.3d(s1,s2,n,H) dev.off()
そんなわけで、outerに投げて動くようにしてみた。逆行列の計算とか自分でやっているため、すごく一般性がないんだけど、二次元plotのみ限定とかだ使える結構速度が出る。
mkr_with_noise<- function(x1,x2){ return(sin(2*pi*x1)+x2+rnorm(1,0,1/4)) } epanechnikov.radial <- function(X1,X2,h11,h12,h21,h22){ X1 <- X1 X2 <- X2 h11 <- h11 h12 <- h12 h21 <- h21 h22 <- h22 return(function(x1,x2){ det <- (h11*h22 - h12*h21) u1 <- 1/det * (h22 * (x1 - X1) - h12 * (x2 - X2)) u2 <- 1/det * (-h21 * (x1 - X1) + h11 * (x2 -X2)) utu <- sqrt(u1^2 + u2^2) ((1-utu) * ifelse(utu <= 1,1,0)) }) } y.exp <- function(X1,X2,Y,h11,h12,h21,h22){ X1 <- X1 X2 <- X2 Y <- Y h11 <- h11 h12 <- h12 h21 <- h21 h22 <- h22 return(function(x1,x2){ det <- (h11*h22 - h12*h21) u1 <- 1/det * (h22 * (x1 - X1) - h12 * (x2 - X2)) u2 <- 1/det * (-h21 * (x1 - X1) + h11 * (x2 -X2)) utu <- sqrt(u1^2 + u2^2) ((1-utu) * ifelse(utu <= 1,1,0)) * Y }) } plot.3d <- function(s1,s2,n,h11,h12,h21,h22){ persp(s1,s2, Reduce("+",Map(function(f){ outer(s1,s2,f) },mapply(function(X1,X2,Y){y.exp(X1,X2,Y,h11,h12,h21,h22)},X1,X2,Y))) / Reduce("+",Map(function(f){ outer(s1,s2,f) },mapply(function(X1,X2){epanechnikov.radial(X1,X2,h11,h12,h21,h22)},X1,X2))), theta=320,phi=20,col=rainbow(50),ticktype="detailed", xlab="x1",ylab="x2",zlab="density") } n <- 30 X1 <- seq(0,1,length=n) X2 <- X1 X <- expand.grid(X1,X2) Y <- apply(X,1,function(X){ x1 <- X[1] x2 <- X[2] mkr_with_noise(x1,x2) }) X1 <- X[,1] X2 <- X[,2] s1 <- seq(0,1,length=n) s2 <- seq(0,1,length=n) setwd("/tmp") png("nadaraya.multi1.png") plot.3d(s1,s2,n,0.2,0,0,0.2) dev.off() png("nadaraya.multi2.png") plot.3d(s1,s2,n,0.1,0,0,0.2) dev.off() png("nadaraya.multi3.png") plot.3d(s1,s2,n,0.1,0.2,0.2,0.1) dev.off() png("nadaraya.multi4.png") plot.3d(s1,s2,n,0.3,0.15,0.1,0.2) dev.off()
