読者です 読者をやめる 読者になる 読者になる

速度を取るか、美しさを取るか

卒研 R

以前多変数ナダラヤワトソン推定量を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()