SVMを作ってみた

この記事は参考にはなりません。コードが汚いのと何も説明を書きません。完全に自己満足用です。

本当はこの記事消したいけど頑張ったからとりあえず置いときます。強くなってまた書き直します。

夏休みに入ったのでsvmの理論をちょいと勉強して実装しました。

###参考資料
[サポートベクトルマシン (機械学習プロフェッショナルシリーズ)](https://www.kspub.co.jp/book/detail/1529069.html)

###SVM本体

収束条件についてはいまいちわからなかった。

```go:SVM.go
package SVM

import (
	"fmt"
	"../matrix"
	"math/rand"
	"math"
	"time"
)

//Cじゃよ
var c = 1.0

//トレーニングデータ（ラベル・目的変数とか呼ばれる方）じゃよ
var label []float64
//トレーニングデータ（特徴量・説明変数とか呼ばれる方）じゃよ
var data [][]float64
//双対変数
var alpha []float64
//KKT条件に違反しているかどうかを判別する
//変数選択に使用するスライス
var ban []float64
//直近で変更した双対変数のインデックス
var index []int
//変更したときに計算したΔβ
var delta_beta float64
//バイアス項
var b float64
//初期化処理じゃよ
func Initdata(C float64,y []float64,x [][]float64){
	c = C
	label = y
	data = x
	alpha = make([]float64,len(label))
	for i:=0;i<len(label);i++{
		//0～Cの間でalphaの初期値とする
		alpha[i] = math.Abs(rand.Float64()) * c
	}
	index = []int{-1,-1}
	rand.Seed(time.Now().UnixNano())
}

//変数選択をする関数
func choice(kernel func([]float64,[]float64) float64){
	//banスライスを更新
	if index[0] >= 0 && index[1] >= 0{
		for i:=0;i<2;i++{
			n := 0.0
			for j:=0;j<len(ban);j++{
				n += label[j]*alpha[j]*kernel(data[index[i]],data[j])
			}
			ban[index[i]] = label[index[i]] - n
		}
		for i:=0;i<len(alpha);i++{
			if i == index[0] || i == index[1]{continue}
			ban[i] -= (kernel(data[i],data[index[0]]) - kernel(data[i],data[index[1]])) * delta_beta
		}
	}
	//0<=α<cの場合と0<α<=cの場合に分けて
	//KKT条件の違反度合いを調べる
	s := 0
	svalue := ban[0]
	for i:=1;i<len(alpha);i++{
		if 0.00001 <= alpha[i] && alpha[i] < c{
			if svalue < ban[i]{
				s = i
				svalue = ban[i]
			}
		}
	}
	t := 0
	if s == 0{t = 1}
	tvalue := ban[t]
	for i:=t+1;i<len(alpha);i++{
		if i == s{continue}
		if 0 <alpha[i] && alpha[i] <= c-0.00001{
			if tvalue > ban[i]{
				t = i
				tvalue = ban[i]
			}
		}
	}
	index[0] = s
	index[1] = t
}

//εを計算する関数
func Epsilon(kernel func([]float64,[]float64) float64) float64 {
	s := index[0];t := index[1]
	x1 := data[s];x2 := data[t]
	y1 := label[s];y2 := label[s]
	a := y1-y2
	for i:=0;i<len(alpha);i++{
		a -= alpha[i]*label[i]*(kernel(data[i],x1) - kernel(data[i],x2))
	}
	b := kernel(x1,x1) - 2*kernel(x1,x2) + kernel(x2,x2)
	return a / b
}

//Δβの範囲を計算する関数
func DeltaBetaRange() (float64,float64) {
	sIndex := index[0];tIndex := index[1]
	if label[sIndex] > 0 && label[tIndex] > 0{
		a := matrix.Max(-alpha[sIndex]*label[sIndex],alpha[tIndex]*label[tIndex]-c)
		b := matrix.Min(c-alpha[sIndex]*label[sIndex],alpha[tIndex]*label[tIndex])
		return a,b
	} else if label[sIndex] > 0 && label[tIndex] < 0{
		a := matrix.Max(-alpha[sIndex]*label[sIndex],alpha[tIndex]*label[tIndex])
		b := matrix.Min(c-alpha[sIndex]*label[sIndex],c+alpha[tIndex]*label[tIndex])
		return a,b
	} else if label[sIndex] < 0 && label[tIndex] > 0 {
		a := matrix.Max(-c-alpha[sIndex]*label[sIndex],alpha[tIndex]*label[tIndex]-c)
		b := matrix.Min(-alpha[sIndex]*label[sIndex],alpha[tIndex]*label[tIndex])
		return a,b
	} else if label[sIndex] < 0 && label[tIndex] < 0 {
		a := matrix.Max(-c-alpha[sIndex]*label[sIndex],alpha[tIndex]*label[tIndex])
		b := matrix.Min(c-alpha[sIndex]*label[sIndex],c+alpha[tIndex]*label[tIndex])
		return a,b
	}
	return 0.0,0.0
}



func SVC(functionName string) ([]float64,float64) {
	f := func([]float64,[]float64) float64 {return 0.0}
	if functionName == "liner"{
		f = func(X_i []float64, X_t []float64) float64 {
			n := 0.0
			for i:=0;i<len(X_i);i++{
				n += X_i[i] * X_t[i]
			}
			return n
		}
	}
	//banの初期化
	ban = make([]float64,len(label))
	for i:=0;i<len(alpha);i++{
		nn := 0.0
		for j:=0;j<len(alpha);j++{
			nn += label[j]*alpha[j]*f(data[i],data[j])
		}
		ban[i] = label[i] - nn
	}

	for i:=0;i<10000;i++ {
		choice(f)
		e := Epsilon(f)
		a,b := DeltaBetaRange()
		//fmt.Println(a,b)
		if a > e {
			delta_beta = a
		} else if e > b{
			delta_beta = b
		} else {
			delta_beta = e
		}
		//fmt.Println(delta_beta)
		//alphaスライスの更新
		alpha[index[0]] = (alpha[index[0]] * label[index[0]] + delta_beta) * label[index[0]]
		alpha[index[1]] = (alpha[index[1]] * label[index[1]] - delta_beta) * label[index[1]]
	}
	var b2,count float64
	count = 0
	w := make([]float64,len(data[0]))
	for i:=0;i<len(data);i++{
		if 0.001 < alpha[i] && alpha[i] < c-0.001{
			count += 1.0
			if count == 1.0{
				for i:=0;i<len(data);i++{
					for j:=0;j<len(data[0]);j++{
						w[j] += data[i][j] * alpha[i] * label[i]
					}
				}
			}
			b2 += label[i] - f(w,data[i])
		}
	}
	fmt.Println(w,b2/count)
	fmt.Println("のん!")
	return w,b2/count
}
```

###SVMを使うためのやつ

```go:main.go
package main

import (
	//"gonum.org/v1/plot"
	//"gonum.org/v1/plot/plotter"
	//"gonum.org/v1/plot/vg"
	"./SVM"
	"math/rand"
	"fmt"
	//"os"
	"os/exec"
	"io"
	"strconv"
	//"time"
)

var x [][]float64
var y []float64

//データの作成
func init()  {
	rand.Seed(0)
	//データの数
	n := 100
	//データの次元数
	d := 2
	//係数行列
	//w := []float64{1,2}
	//バイアス
	//b := 3.0
	//クラス
	class := [][]float64{{4,4},{3,3}}
	x = make([][]float64,n)
	y = make([]float64,n)
	for i:=0;i<n;i++{
		x[i] = make([]float64,d)
		for j:=0;j<2;j++{
			k := 0
			if n/2 < i{k = 1}
			x[i][j] = rand.Float64() + class[k][j]
		}
		if n/2 < i{
			y[i] = 1
		} else {
			y[i] = -1
		}
		if rand.Intn(101) < 10{
			if y[i] == -1{
				y[i] = 1
			} else {
				y[i] = -1
			}
		}
	}
}

func main()  {
	SVM.Initdata(0.1,y,x)
	w,b := SVM.SVC("liner")
	//fmt.Println(w)
	//fmt.Println(b)

	cmd := exec.Command("python","C:/pythonfile/machen_larning/スクラッチ実装/SVM/plot.py")
	stdin,_ := cmd.StdinPipe()
	io.WriteString(stdin,strconv.Itoa(len(x)) + "\n")
	for i:=0;i<len(x);i++{
		io.WriteString(stdin,strconv.FormatFloat(x[i][0], 'f', 15, 64) + " " + strconv.FormatFloat(x[i][1], 'f', 15, 64) + "\n")
	}
	for i:=0;i<len(y);i++{
		io.WriteString(stdin,strconv.Itoa(int(y[i])) + " ")
	}
	io.WriteString(stdin,"\n")
	io.WriteString(stdin, strconv.FormatFloat(b, 'f', 15, 64) + "\n")
	for i:=0;i<len(w);i++{
		io.WriteString(stdin, strconv.FormatFloat(w[i], 'f', 15, 64) + " ")
	}
	io.WriteString(stdin,"\n")
	stdin.Close()
	out, _ := cmd.CombinedOutput()
	fmt.Println(string(out))
}
```

###呼び出すpythonのプログラム

グラフ書く用です。

```python:plot.py
from matplotlib import pyplot as plt

n = int(input())
x1 = [0 for i in range(n)];x2 = [0 for i in range(n)]
for i in range(n):
    x1[i],x2[i] = list(map(float,input().split()))
y = [int(i) for i in input().split()]
#マイナスの時とプラスの時で分ける
x1min = min(x1);x1max = max(x1)
xm = [[],[]];xp = [[],[]]
for i in range(len(y)):
    if y[i] < 0:
        xm[0].append(x1[i])
        xm[1].append(x2[i])
    else:
        xp[0].append(x1[i])
        xp[1].append(x2[i])


b = float(input())
w = [float(i) for i in input().split()]
line = [-b/i for i in w]
x = [i for i in range(int(x1min-1),int(x1max+1)+2)]
y = [-(b+i*w[0])/w[1] for i in x]

fig = plt.figure()
ax = fig.add_subplot(1,1,1)
ax.scatter(xm[0],xm[1], c='red')
ax.scatter(xp[0],xp[1], c='blue')
ax.plot(x,y)
ax.set_xlim([x1min-1, x1max+1])
plt.show()
```


![Figure_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/278116/cf639f9f-f399-4938-78e3-2f24dde06428.png)

#さいごに

疲れました。