第8章　サポートベクトルマシン（問題75～87）

82

library(quadprog)

# 関数の定義
svm.1 <- function(X, y, C) {
  eps <- 0.0001
  n <- nrow(X)
  Dmat <- matrix(nrow = n, ncol = n)
  for (i in 1:n)
    for (j in 1:n)
      Dmat[i, j] <- sum(X[i, ] * X[j, ]) * y[i] * y[j]
  Dmat <- Dmat + eps * diag(n)
  dvec <- ## 空欄(1) ##
  Amat <- matrix(nrow = (2 * n + 1), ncol = n)
  Amat[1, ] <- y
  Amat[2:(n + 1), 1:n] <- -diag(n)
  Amat[(n + 2):(2 * n + 1), 1:n] <- ## 空欄(2) ##
  Amat <- t(Amat)
  bvec <- ## 空欄(3) ##
  meq <- 1
  alpha <- solve.QP(Dmat, dvec, Amat, bvec = bvec, meq = 1)$solution
  beta <- drop((alpha * y) %*% X)
  index <- (1:n)[eps < alpha & alpha < C - eps]
  beta.0 <- mean(y[index] - X[index, ] %*% beta)
  return(list(beta = beta, beta.0 = beta.0[1]))
}

# svm.1を用いた実行
a <- rnorm(1)
b <- rnorm(1)
n <- 100
X <- matrix(rnorm(n * 2), ncol = 2, nrow =n)
y <- sign(a * X[, 1] + b * X[, 2] + 0.3 * rnorm(n))
plot(-3:3, -3:3, xlab = "X[, 1]", ylab = "X[, 2]", type = "n")
for (i in 1:n) {
  if (y[i] == 1) {
    points(X[i, 1], X[i, 2], col = "red")
  } else {
    points(X[i, 1], X[i, 2], col = "blue")
  }
}
qq <- svm.1(X, y, 10)
abline(-qq$beta.0 / qq$beta[2], -qq$beta[1] / qq$beta[2])

85

# 関数の定義
plot.kernel <- function(K, lty) {
  qq <- svm.2(X, y, 1, K)
  alpha <- qq$alpha
  beta.0 <- qq$beta.0

  f <- function(u, v) {
    x <- c(u, v)
    S <- beta.0
    for (i in 1:n)
      S <- S + ## 空欄 ##
    return(S)
  }
  
  u <- seq(-2, 2, .1)
  v <- seq(-2, 2, .1)
  w <- array(dim = c(41, 41))
  for (i in 1:41)
    for (j in 1:41)
      w[i, j] <- f(u[i], v[j])
  contour(u, v, w, level = 0, add = TRUE, lty = lty)
}

# 実行
a <- 3
b <- -1
n <- 200
X <- matrix(rnorm(n * 2), ncol = 2, nrow = n)
y <- sign(a * X[, 1] + b * X[, 2] ^ 2 + 0.3 * rnorm(n))
plot(-3:3, -3:3, xlab = "X[, 1]", ylab = "X[, 2]", type = "n")
for (i in 1:n) {
  if (y[i] == 1) {
    points(X[i, 1], X[i, 2], col = "red")
  } else {
    points(X[i, 1], X[i, 2], col = "blue")
  }
}
plot.kernel(K.linear, 1)
plot.kernel(K.poly, 2)

86

library(e1071)

x <- matrix(rnorm(200 * 2), ncol = 2)
x[1:100, ] <- x[1:100, ] + 2
x[101:150, ] <- x[101:150, ] - 2
y <- c(rep(1, 150), rep(2, 50))
dat <- data.frame(x = x, y = as.factor(y))
train <- sample(200, 100)
svmfit <- svm(y ~ ., data = dat[train, ], kernel = "radial",
              gamma = 1, cost = 1)
plot(svmfit, dat[train, ])

tune.out <- tune(svm, y ~ ., data = dat[train, ], kernel = "radial",
                 ranges = list(cost = c(0.1, 1, 10, 100, 1000),
                               gamma = c(0.5, 1, 2, 3, 4)))
summary (tune.out)

87

library(e1071)

df <- iris
x <- df[, 1:4]
y <- as.factor(df[, 5])
m <- 30
train <- sample(1:150, 150 - m)
x.train <- x[train, ]
y.train <- y[train]
dat <- data.frame(## 空欄(1) ##, y.train)
svmfit <- svm(## 空欄(2) ## ~ ., data = dat, kernel = "radial",
              cost = 10, gamma = 1)
x.test <- x[-train, ]
y.test <- y[-train]
table(y.test, predict(svmfit, x.test))