PRML 7.1.1に記載の、一部のデータの誤分類を許すようなサポートベクトルマシンを構成し、いくつかのCの値毎に、推定値、および各サンプルのラグランジュ乗数aを図示します。なお、カーネル関数は用いておらず、線形カーネル相当となります。二次計画法には、kernlabパッケージのipopを用いています。
library(MASS)
library(kernlab)
frame()
set.seed(0)
par(mfrow=c(2, 2))
par(mar=c(3, 3, 1, 0.1))
par(mgp=c(2, 1, 0))
xrange <- c(-4, 8)
yrange <- c(-8, 4)
x1 <- seq(xrange[1], xrange[2], .1)
x2 <- seq(yrange[1], yrange[2], .1)
x <- mvrnorm(10, c(0,1), matrix(c(2, 1, 1, 2), 2))
d1 <- as.data.frame(cbind(x, 1, 1))
x <- mvrnorm(10, c(2,0), matrix(c(2, 1, 1, 2), 2))
d1 <- rbind(d1, as.data.frame(cbind(x, 2, -1)))
x <- mvrnorm(5, c(7,-7), matrix(c(.2, 0, 0, .2), 2))
d2 <- rbind(d1, as.data.frame(cbind(x, 2, -1)))
d3 <- as.data.frame(rbind(cbind(1, 1, 1, 1), cbind(0, 1, 2, -1)))
names(d1) <- c("x1","x2","class","t")
names(d2) <- c("x1","x2","class","t")
names(d3) <- c("x1","x2","class","t")
doplot <- function(d, C) {
N <- nrow(d)
# quadratic programming
cc <- matrix(rep(-1, N))
hmat <- cbind(d$t * d$x1, d$t * d$x2)
hmat <- hmat %*% t(hmat)
ll <- matrix(rep(0, N))
uu <- matrix(rep(C, N))
bb <- 0
rr <- 0
amat <- matrix(d$t, nrow=1, ncol=N)
sv <- ipop(cc, hmat, amat, bb, ll, uu, rr)
a <- primal(sv)
# non-support vectors (a=0)
a[a < 1.0e-5] <- 0
# support vectors with penalty (a=C)
a[a > C - 1.0e-5] <- C
cat("a");print(a)
# derive b
bs <- c()
for (n in 1:N) {
if (a[n] > 0 && a[n] < C) {
bs <- c(bs, d$t[n] - rowSums(a * d$t * (cbind(d$x1[n], d$x2[n]) %*% t(cbind(d$x1, d$x2)))))
}
}
b <- mean(bs)
# plot estimates
estimate <- function(x1, x2) {
result <- 0
for (n in 1:N) {
if (a[n] != 0) {
result <- result + a[n] * d$t[n] * (cbind(x1, x2) %*% t(cbind(d$x1[n], d$x2[n])))
}
}
result <- result + b
}
y1 <- outer(x1, x2, estimate)
plot(d$x1, d$x2, col=d$class, type="n", xlim=xrange, ylim=yrange)
image(x1, x2, y1, xlim=xrange, ylim=yrange, col=terrain.colors(32), add=T)
contour(x1, x2, y1, xlim=xrange, ylim=yrange, levels=c(-1, 0, 1), add=T)
par(new=T)
plot(d$x1, d$x2, pch=ifelse(a == 0, 1, ifelse(a == C, 17, 16)), col=d$class, xlim=xrange, ylim=yrange, cex=1.3)
legend("bottomleft", legend=c("a=0","0<a<C","a=C"), pch=c(1, 16, 17), cex=1.3, y.intersp=1, bg="white")
title(paste("C=", C))
}
doplot(d2, 0.5)
doplot(d2, 1)
doplot(d2, 2)
doplot(d2, 4)