入門 機械学習による異常検知

3.3.4 節（カーネル密度推定）

実行例3.5

library(car)
library(KernSmooth)
x <- Davis[,c('weight','height')]
h <- c(dpik(x$weight), dpik(x$height))
est <- bkde2D(x,bandwidth=h, gridsize=c(10^3,10^3))
d <- list(x=est$x1,y=est$x2, z=est$fhat)
image(d,col=terrain.colors(7),xlab="weight",ylab="height",xlim=c(35,110),ylim=c(145,200))
contour(d,add=T)

カーネル行列の計算

この部分が本には書かれておりませんでした。色々編集を行っているうちに抜けてしまったものです。失礼いたしました。

# kernel matrix (missing - sorry!)
n <- nrow(x); K <- matrix(-1,n,n)
prefac <- (2*pi*h)^(-0.5)
for(nn in 1:n){
  xnn <- x[nn,]
  x_x1 <- x[,1] - as.numeric(xnn[1])
  K1 <- prefac[1]*exp(-0.5*(1/h[1])^2*(x_x1*x_x1))
  x_x2 <- x[,2] - as.numeric(xnn[2])
  K2 <- prefac[2]*exp(-0.5*(1/h[2])^2*(x_x2*x_x2))
  K[,nn] <- K1 * K2
}

実行例 3.6

# Code 3.6: anomaly scores
aa <- colSums(K)-diag(K)
lowerLim <- 10^(-20); aa[(aa<lowerLim)] <- lowerLim
a <- (-1)*log(aa/(n-1))
plot(a,xlab="sample ID",ylab="anomaly score")