library(MASS)
n <- 100
c <- vector('list', 7)
set.seed(777)
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')
}

# 単純な分類用データ(クラスターサイズ:2)
d2 <- data.frame(c(rep(1, n), rep(0, n)), rbind(c[[1]], c[[5]]))
colnames(d2) <- c('blue', 'x', 'y')
write.csv(d2, file = 'data_svm_cluster2.csv', row.names = F)

# 複雑な分類用データ(クラスターサイズ:7)
d7 <- data.frame(c(rep(1, n*3), rep(0, n*4)),
                  rbind(c[[1]], c[[2]], c[[3]], c[[4]], c[[5]], c[[6]], c[[7]]))
colnames(d7) <- c('blue', 'x', 'y')
write.csv(d7, file = 'data_svm_cluster7.csv', row.names = F)
# カラーパレット
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)) # 灰
# グラフ作成関数(後で再利用するため関数化)
draw.fig <- function(d2)
{
  # データ抽出
  d.red  <- d2[d2$blue == 0, ] # 赤クラスデータ
  d.blue <- d2[d2$blue == 1, ] # 青クラスデータ

  # 図枠
  matplot (NA, type = 'n',
           xlim = c(-10, 15), ylim = c(-10, 20),
           xlab = 'x', ylab = 'y')

  grid() # 格子線 

  # 描画
  matlines(x = d.red$x,  y = d.red$y,  type = 'p', pch = 1, col = COL[2])
  matlines(x = d.blue$x, y = d.blue$y, type = 'p', pch = 1, col = COL[1])

  # 凡例
  legend('topright', col = COL[1:2], pch = c(1, 1), bg = 'white',
        legend = c('赤', '青'))
}

draw.fig(d2)

# テストデータ読込
d2 <- read.csv(file = 'https://stats.dip.jp/01_ds/data/data_svm_cluster2.csv')

# カーネル関数
KERNEL <- c('linear',     # 線形
            'polynomial', # 多項式
            'sigmoid',    # シグモイド
            'radial')     # ガウス

k <- 1 # カーネル選択番号

# 交差検証法によるパラメータ選択
library(e1071)
cv <- tune('svm', as.factor(blue) ~ ., data = d2,
           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:
##  cost
##     2
## 
## - best performance: 0.01
dx <- 0.2
dy <- 0.2

# 格子点データを作成
dgrid <- expand.grid(x = seq(-25, 25, dx),
                     y = seq(-25, 25, dy))

# 格子点を予測
pred <- predict(cv$best.model, newdata = dgrid)

draw.fig(d2)

# サポートベクター
sv <- d2[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])
}

# 予測値が1(青)の周りを正方形で灰色塗り(ポリゴンで埋める)
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('赤(0)', '青(1)', 'サポートベクター', '青(1)と分類する範囲'))

KERNEL # カーネル関数の種類
## [1] "linear"     "polynomial" "sigmoid"    "radial"
d2 <- read.csv(file = 'https://stats.dip.jp/01_ds/data/data_svm_cluster2.csv')

# カーネル関数
KERNEL <- c('linear',     # 線形
            'polynomial', # 多項式
            'sigmoid',    # シグモイド
            'radial')     # ガウス

k <- 2 # カーネル選択番号

# 交差検証法によるパラメータ選択
library(e1071)
cv <- tune('svm', as.factor(blue) ~ ., data = d2,
           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:
##  coef0 cost
##  0.125    8
## 
## - best performance: 0
dx <- 0.2
dy <- 0.2

# 格子点データを作成
dgrid <- expand.grid(x = seq(-25, 25, dx),
                     y = seq(-25, 25, dy))

# 格子点を予測
pred <- predict(cv$best.model, newdata = dgrid)

draw.fig(d2)

# サポートベクター
sv <- d2[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])
}

# 予測値が1(青)の周りを正方形で灰色塗り(ポリゴンで埋める)
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('赤(0)', '青(1)', 'サポートベクター', '青(1)と分類する範囲'))

d2 <- read.csv(file = 'https://stats.dip.jp/01_ds/data/data_svm_cluster2.csv')

# カーネル関数
KERNEL <- c('linear',     # 線形
            'polynomial', # 多項式
            'sigmoid',    # シグモイド
            'radial')     # ガウス

k <- 3 # カーネル選択番号

# 交差検証法によるパラメータ選択
library(e1071)
cv <- tune('svm', as.factor(blue) ~ ., data = d2,
           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:
##  coef0 cost
##      4    2
## 
## - best performance: 0.01
dx <- 0.2
dy <- 0.2

# 格子点データを作成
dgrid <- expand.grid(x = seq(-25, 25, dx),
                     y = seq(-25, 25, dy))

# 格子点を予測
pred <- predict(cv$best.model, newdata = dgrid)

draw.fig(d2)

# サポートベクター
sv <- d2[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])
}

# 予測値が1(青)の周りを正方形で灰色塗り(ポリゴンで埋める)
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('赤(0)', '青(1)', 'サポートベクター', '青(1)と分類する範囲'))

d2 <- read.csv(file = 'https://stats.dip.jp/01_ds/data/data_svm_cluster2.csv')

# カーネル関数
KERNEL <- c('linear',     # 線形
            'polynomial', # 多項式
            'sigmoid',    # シグモイド
            'radial')     # ガウス

k <- 4 # カーネル選択番号

# 交差検証法によるパラメータ選択
library(e1071)
cv <- tune('svm', as.factor(blue) ~ ., data = d2,
           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
##      8    1
## 
## - best performance: 0.005
dx <- 0.2
dy <- 0.2

# 格子点データを作成
dgrid <- expand.grid(x = seq(-25, 25, dx),
                     y = seq(-25, 25, dy))

# 格子点を予測
pred <- predict(cv$best.model, newdata = dgrid)

draw.fig(d2)

# サポートベクター
sv <- d2[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])
}

# 予測値が1(青)の周りを正方形で灰色塗り(ポリゴンで埋める)
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('赤(0)', '青(1)', 'サポートベクター', '青(1)と分類する範囲'))

d7 <- 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() # 格子線

# データ抽出
d7.red  <- d7[d7$blue == 0, ] # 赤データ
d7.blue <- d7[d7$blue == 1, ] # 青データ

# 描画
matlines(x = d7.red$x,  y = d7.red$y,  type = 'p', pch = 1, col = COL[1])
matlines(x = d7.blue$x, y = d7.blue$y, type = 'p', pch = 1, col = COL[2])

# 凡例
legend('topright', col = COL[1:2], pch = c(1, 1), bg = 'white',
      legend = c('赤', '青'))

# グラフ作成関数(後で再利用するため関数化)
draw.fig <- function(d7)
{
  # データ抽出
  d.red  <- d7[d7$blue == 0, ] # 赤クラスデータ
  d.blue <- d7[d7$blue == 1, ] # 青クラスデータ

  # 図枠
  matplot (NA, type = 'n',
           xlim = c(-10, 15), ylim = c(-10, 20),
           xlab = 'x', ylab = 'y')

  grid() # 格子線 

  # 描画
  matlines(x = d.red$x,  y = d.red$y,  type = 'p', pch = 1, col = COL[2])
  matlines(x = d.blue$x, y = d.blue$y, type = 'p', pch = 1, col = COL[1])

  # 凡例
  legend('topright', col = COL[1:2], pch = c(1, 1), bg = 'white',
        legend = c('赤', '青'))
}

draw.fig(d7)

# カーネル関数
KERNEL <- c('linear',     # 線形
            'polynomial', # 多項式
            'sigmoid',    # シグモイド
            'radial')     # ガウス

k <- 4 # カーネル選択番号

# 交差検証法によるパラメータ選択
library(e1071)
cv <- tune('svm', as.factor(blue) ~ ., data = d7,
           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
##      2    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(d7)

# サポートベクター
sv <- d7[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])
}

# 予測値が1(青)の周りを正方形で灰色塗り(ポリゴンで埋める)
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('赤(0)', '青(1)', 'サポートベクター', '青(1)と分類する範囲'))