library(MASS)
n <- 100
c <- vector('list', 7)
c[[1]] <- mvrnorm(n, mu = c( 0, 0), Sigma = rbind(c(2.0, 0.0), c( 0.0, 2.0)))
c[[2]] <- mvrnorm(n, mu = c( 0, 10), Sigma = rbind(c(2.0, -0.8), c(-0.8, 2.0)))
c[[3]] <- mvrnorm(n, mu = c(10, 0), Sigma = rbind(c(2.0, -0.8), c(-0.8, 2.0)))
c[[4]] <- mvrnorm(n, mu = c(-5, -5), Sigma = rbind(c(2.0, 0.8), c( 0.8, 2.0)))
c[[5]] <- mvrnorm(n, mu = c( 5, 5), Sigma = rbind(c(2.0, 0.8), c( 0.8, 2.0)))
c[[6]] <- mvrnorm(n, mu = c(-5, 5), Sigma = rbind(c(2.0, -0.8), c(-0.8, 2.0)))
c[[7]] <- mvrnorm(n, mu = c( 5, -5), Sigma = rbind(c(2.0, -0.8), c(-0.8, 2.0)))
for (i in seq_along(c))
{
c[[i]] <- as.data.frame(c[[i]])
colnames(c[[i]]) <- c('x', 'y')
}
COL <- c(rgb(255, 0, 0, 105, max = 255), # 赤
rgb( 0, 0, 255, 105, max = 255), # 青
rgb( 0, 155, 0, 105, max = 255), # 緑
rgb(100, 100, 100, 20, max = 255)) # 灰
# 単純な分類用データ
d <- data.frame(c(rep(1, n), rep(0, n)), rbind(c[[1]], c[[5]]))
colnames(d) <- c('blue', 'x', 'y')
draw.fig <- function()
{
# 図枠
matplot (NA, type = 'n',
xlim = c(-10, 15), ylim = c(-10, 20),
xlab = 'x', ylab = 'y')
grid() # 格子線
# 描画
matlines(x = c[[1]]$x, y = c[[1]]$y, type = 'p', pch = 1, col = COL[2])
matlines(x = c[[2]]$x, y = c[[2]]$y, type = 'p', pch = 1, col = COL[2])
matlines(x = c[[3]]$x, y = c[[3]]$y, type = 'p', pch = 1, col = COL[2])
matlines(x = c[[4]]$x, y = c[[4]]$y, type = 'p', pch = 1, col = COL[1])
matlines(x = c[[5]]$x, y = c[[5]]$y, type = 'p', pch = 1, col = COL[1])
matlines(x = c[[6]]$x, y = c[[6]]$y, type = 'p', pch = 1, col = COL[1])
matlines(x = c[[7]]$x, y = c[[7]]$y, type = 'p', pch = 1, col = COL[1])
}
#cairo_pdf('data_svm.pdf') # PDF画像作成（ここから）
draw.fig()
# 凡例
legend('topright', col = COL[1:2], pch = c(1, 1), bg = 'white',
legend = c('赤', '青'))
KERNEL <- c('linear', # 線形
'polynomial', # 多項式
'sigmoid', # シグモイド
'radial') # ガウス
k <- 4
d0 <- read.csv(file = 'https://stats.dip.jp/01_ds/data/data_svm_cluster7.csv')
# 図枠
matplot(NA, type = 'n',
xlim = c(-10, 15), ylim = c(-10, 20),
xlab = 'x', ylab = 'y')
grid() # 格子線
# データ抽出
d0.red <- d0[d0$blue == 0, ] # 赤データ
d0.blue <- d0[d0$blue == 1, ] # 青データ
# 描画
matlines(x = d0.red$x, y = d0.red$y, type = 'p', pch = 1, col = COL[1])
matlines(x = d0.blue$x, y = d0.blue$y, type = 'p', pch = 1, col = COL[2])
# 凡例
legend('topright', col = COL[1:2], pch = c(1, 1), bg = 'white',
legend = c('赤', '青'))
library(e1071)
cv <- tune('svm', as.factor(blue) ~ ., data = d0,
kernel = KERNEL[k], type = 'C-classification',
ranges = list(gamma = 2^(-4:4), # radialなどの非線形カーネルを使うとき調整
#epsilon = seq(0, 1, 0.1), # SVRの不感帯の調整
#coef0 = 2^(-4:4), # polynomialかsigmoidのとき調整（c0）
cost = 2^(-4:4))) # コスト係数（小さいほど分類誤りを許容）
# ベストパラメータ表示
cv
##
## Parameter tuning of 'svm':
##
## - sampling method: 10-fold cross validation
##
## - best parameters:
## gamma cost
## 1 4
##
## - best performance: 0.01142857
dx <- 0.1
dy <- 0.1
# 格子点データを作成
dgrid <- expand.grid(x = seq(-25, 25, dx),
y = seq(-25, 25, dy))
# 格子点を予測
pred <- predict(cv$best.model, newdata = dgrid)
draw.fig()
# サポートベクター
sv <- d0[cv$best.model$index, -1]
matpoints(x = sv[, 1], y = sv[, 2], pch = 16, cex = 0.5, col = 1)
dgrid.blue <- dgrid[pred == 1, ]
#matpoints(x = dgrid.blue$x, y = dgrid.blue$y, pch = 15, cex = 1.1, col = COL[4])
# 灰色塗り（ポリゴン）関数
fill.cell <- function(x, y)
{
xline <- c(x - dx/2, x + dx/2)
ylow <- c(y - dy/2, y - dy/2)
yupp <- c(y + dy/2, y + dy/2)
polygon(c(xline, rev(xline)), c(ylow, yupp), border = F, col = COL[4])
}
# 予測値が１（青）の周りを正方形で灰色塗り（ポリゴンで埋める）
for (i in 1:nrow(dgrid))
{
if (pred[i] == 1) fill.cell(dgrid$x[i], dgrid$y[i])
}
# 主タイトル
title(paste0('SVM（カーネル：', KERNEL[k], '）による分類'))
# 凡例
legend('topright', col = c(COL[1:2], 1, NA), pch = c(1, 1, 16, NA),
fill = c(NA, NA, NA, COL[4]), border = F, bg = 'white',
legend = c('赤（０）', '青（１）', 'サポートベクター', '青（１）と分類する範囲'))
