Business Analytics

30 Oct 2018

Support Vector Machine

고려대학교 Business Analytics 수업의 과제를 수행하기 위해 제작하였습니다 -
해당 포스트는 전적으로 고려대학교 강필성 교수님의 Business Analytics 강의에 바탕을 둡니다.-

SVM 은 서로 다른 클래스의 값들을 분류(classify)하기 위해 사용하는 알고리즘입니다. 아래의 그림에서와 같이 실선을 통해 별표와 원을 구분할 수 있을 것입니다.

이 때 실선과 점선 사이의 거리(Gap)을 ‘마진(margin)’이라고 합니다. ‘마진(margin)’은 SVM 에서 가장 중요한 개념 중 하나입니다. 이에 대해서는 차차 알아가도록 하겠습니다.

아래의 그림을 보면 분류모델 B1, B2 두 가지가 있는 것을 알 수 있습니다. 이 때 어떤 분류선이 더 좋은 선일까요? 알고리즘마다 다른 기준이 적용됩니다. 그러나 SVM에 국한해서 봤을 때, 더 좋은 분류선은 마진(margin)이 더 큰 B1 입니다.

모델이 똑같은 정확도를 가진다고 했을 때, 마진(margin)이 클 수록 VC-dimension은 작게되고, 이에 따라 모델의 구조적 위험성도 줄어들게 되기 때문입니다.
따라서 앞으로의 목적은 마진(margin)을 최대화하는 분류선을 찾는 것이 될 것입니다.

Margin

아래의 그림은 데이터를 분류하는 과정을 3차원으로 나타낸 그림입니다. 해당 plus-plane 보다 위에 있는 데이터는 plus 로 분류되고, minus-plane 보다 아래에 있는 데이터는 minus로 분류됩니다.
따라서 plus 에 해당하는 데이터는 $w^tx + b \ge 1$ 를 만족하고, minus에 해당하는 데이터는 $w^tx + b \le -1$를 만족합니다.

수식에 담긴 의미를 조금 더 해석해보겠습니다. $w$ 는 plane 들을 직교하며 지나는 normale vector 이고, x는 데이터벡터입니다. $w^t x$는 normale vector와 데이터 간의 내적으로, 그 값이 1보다 크면 plus value 가 되고, 값이 -1보다 작으면 minus value가 되는 것입니다.

SVM 은 내적이라는 단순한 방법을 통해 데이터를 분류한다는 것을 알 수 있는 대목입니다.

Marging 은 두 분류경계간의 거리이므로, plus plane과 minus plane 간의 거리가 곧 Margin 입니다. 이제 $w$를 통해 Margin 을 찾을 차례입니다.

위의 그림을 통해 $x_+$
는 $x_-$ 가 ‘normal vector’인 $w$ 방향으로 p(미지수)만큼 이동했을 때의 값인 것을 알 수 있습니다. 이를 수식화하겠습니다.

$\overset{\to}{x}_{+} = \overset{\to}{x}_{-} + p*w \qquad (1)$

$x^{+}$ 와 $x^{-}$ 에 대한 관계식이 나왔으니 이제 대입만 하면 됩니다.

$w^{T}x + b = 1 \qquad (2)$ $w^{T}x^{+} + b = 1 \qquad (3)$

(1)을 (3)에 대입하면

$w^{T}\left( x^{-} + pw \right) + b = 1$ $w^{T}x^{-} + b + pw = 1$ $-1 + pw^{T}w = 1$ $pw{T}w = 2$ $p = \frac {2} {w{T}w}$

위의 식을 차례대로 전개할 수 있게 됩니다.
p는 normal vector 방향으로 이동하는 값이므로 plane 간의 거리를 뜻합니다.
즉 p값이 margin 값이 되고, 이는 $ \frac {2} {w{T}w}$ 이 되는 것입니다

$Margin = p = \frac {2} {w{T}w}$

앞서 언급했듯이 Margin 이 클 수록 모델이 강건해지므로 Margin 을 최대화해야합니다. Margin 에 대한 최적화는 뒤에서 알아보겠습니다.

Optimization (Maximize Margin)

마진을 최대화하기 위해 $\frac {2} {w{T}w}$ 을 maximize 하는 것은 계산이 불편할 것 같습니다.
편리한 계산을 위해 목적식에 역수를 취하는 변형을 주도록 합니다.
최대화하는 문제에 역수를 취했으므로, 그 문제는 목적이 최소화를 하는 것이라는 유추를 할 수 있습니다.

$max \frac {2} {w{T}w} \Rightarrow min \frac {1} {2} \lVert {w} \rVert^{2}$

제약식은 마진을 이루는 분류경계면이 됩니다.
하나의 식에 두 가지를 넣기 위해 $y_i$ 로 트릭을 줍니다.

$\left( w^{t}x_{i} + b \right) \; \ge 1 \quad , \left( w^{t}x_i +b \right) \; \le -1$

y에 대해 아래와 같이 정의한 후

$y = \begin{cases} 1, & \mbox{for } \mbox{ +} \\ {-1}, & \mbox{for } \mbox{ -} \end{cases}$

분류경계면에 대해 적용해주면 아래와 같은 제약식을 얻을 수 있습니다.

$y_i \left( w^{t}x_i +b \right) \ge 1$

헷갈리는 분을 위해 정리합니다.

$min \frac {1} {2} \lVert {w} \rVert^{2}$ $s.t. \quad y_i \left( w^{t}x_i +b \right) \ge 1$

목적식과 제약식을 구했습니다. 하지만 이 두 개의 식으로 어떻게 해를 구할 수 있을까요? 문제를 더 명료화시키기 위해서 Largrangian Multiplier 를 사용합니다 Largrangian Multiplier 는 제약식에 미지수(여기서는 알파사용)를 곱하고 목적식에 집어 넣어 최적화를 쉽게 해주는 방법입니다. 이 방법 덕분에 우리는 제약식에 신경쓰지 않고 최적화 문제를 풀 수 있게 되었습니다.

목적식과 제약식을 합쳐 라그랑지안 방정식을 만들면 다음과 같은 결과가 나옵니다.

${\min{ L_{p}(w,b,{ \alpha }_{ i }) } } =\frac { 1 }{ 2 } { \left| w \right| }^{ 2 }-\sum { i=1 }^{ n }{ { \alpha }_{ i }({ y }_{ i } ({ w }^{ T }{ x }_{ i }+b)-1) }$

라그랑지안 문제에도 제약은 있습니다.

${ \alpha }_{ i }\ge 0,\quad i=1,...,n$

이런 제약식은 감당 가능합니다.

Dual Problem 으로 변환

라그랑지안만으로도 해당 문제를 푸는 것은 힘듭니다. Dual problem 으로 변환해서 우리가 가진 Lagrangian problem(minimize) 에 대한 lower bound를 생성하고, lower bound 를 maximize 하면 Lagrangian 문제의 최솟값을 구할 수 있습니다.

이를 위해서 $L_p$ 에 대해서 여러 미지수로 각각 편미분한 값을 만들어냅니다. 다들 아시겠지만 편미분한 식이 0이 되면, 그 지점에서 최솟값(lower bound 생성위해)을 구할 수 있습니다. 각 미분값은 다음과 같습니다.

$\begin{align*} \frac { \partial L(w,b,{ \alpha }{ i }) }{ \partial w } =0\quad &\rightarrow \quad w=\sum { i=1 }^{ n }{ { \alpha }{ i }{ y }{ i }{ x }{ i } } \ \frac { \partial L(w,b,{ \alpha }{ i }) }{ \partial b } =0\quad &\rightarrow \quad \sum { i=1 }^{ n }{ { \alpha }{ i }{ y }_{ i } } =0 \end{align*}$

이 미분값들을 $L_p \left(w,b,\alpha_{i} \right)$ 에 대입할 수 있습니다. 대입하여 전개하는 과정은 생략하도록 하겠습니다. 그 결과로 생성되는 Dual Problem 은 다음과 같습니다.

$\max { { L }_{ D }({ \alpha }_{ i }) } =\sum_{ i=1 }^{ n }{ { \alpha }_{ i } } -\frac { 1 }{ 2 } \sum_{ i=1 }^{ n }{ \sum_{ j=1 }^{ n }{ { \alpha }_{ i }{ { \alpha }_{ j }y }_{ i }{ y }_{ j }{ x }_{ i }^{ T }{ x }_{ j } } }$ $s.t. \quad \sum_{ i=1 }^{ n }{ { \alpha }_{ i }{ y }_{ i } } = 0 { \alpha }_{ i }\ge 0,\quad i=1,...,n$

마지막 식이 결국 $\alpha$에 대한 maximization 문제로 정리가 되었습니다. 이 문제를 풀어서 $\alpha$를 구하면, $w= \sum {i=1}^{n}{ \alpha_{i}y_{i}x_{i} }$ 에 대입하여 $w$ 값을 구할 수 있게 됩니다. 뒤이어 도출한 값들과 $y_i \left( {w^{t}x_i +b} \right)-1$ 을 통해서 b역시 구할 수 있게 됩니다.

결과

KKT condition에 의해

${\alpha}_{ i } \left( { y }_{ i }({ w }^{ T }{ x }_{ i }+b)-1 \right) = 0 \quad (4)$

을 도출할 수 있습니다.

(4)를 통해서 $\alpha_{i}$ 또는 $y_i(w^Tx_i+b-1)$ 가 무조건 0이어야 한다는 것을 알 수 있습니다. $\alpha$ = 0이 아니면 $y_i(w^Tx_i+b-1)$ 가 반드시 0입니다.
따라서 $\overset{\rightarrow}{x}$ 는 분류경계면(plane) 위에 있는 벡터가 됩니다.
$\alpha$ = 0 이면 $\overset{\rightarrow}{x}$ 는 분류경계면과 떨어져 있는 벡터가 됩니다

이렇듯 $\alpha$ 가 마진에 대한 영향을 미치므로, $\alpha$를 support vector 라고 합니다.

파이썬 코드를 통해서 위의 내용을 좀 더 자세히 알아보도록 하겠습니다.

import numpy as np
from numpy import linalg
import cvxopt
import cvxopt.solvers
import pylab as pl

사용된 라이브러리는 다음과 같습니다. cvxopt는 최적화를 파이썬으로 구현할 수 있게 도와주는 라이브러리입니다. pylab은 python 통해서 시각화를 할 수 있게하는 라이브러리입니다. 앞서서 최적화 문제를 Lagrangian multiplier, Dual 화를 통해서 풀기 쉽게 유도했기 때문에 그 과정을 코드에서 반복할 필요는 없습니다. 따라서 Dual 대한 수식을 바로 구현하여 $\alpha$ 값을 구합니다. fit 으로 해당하는 과정을 묶어줍니다.

def fit(self, X, y):

    n_samples, n_features = X.shape

    # Gram matrix
    K = np.zeros((n_samples, n_samples))
    for i in range(n_samples):
        for j in range(n_samples):
            K[i,j] = self.kernel(X[i], X[j])

    P = cvxopt.matrix(np.outer(y,y) * K)
    q = cvxopt.matrix(np.ones(n_samples) * -1)
    A = cvxopt.matrix(y, (1,n_samples))
    b = cvxopt.matrix(np.zeros(1))

    if self.C is None:
        G = cvxopt.matrix(np.diag(np.ones(n_samples) * -1))
        h = cvxopt.matrix(np.zeros(n_samples))
    else:
        tmp1 = np.diag(np.ones(n_samples) * -1)
        tmp2 = np.identity(n_samples)
        G = cvxopt.matrix(np.vstack((tmp1, tmp2)))

        tmp1 = np.zeros(n_samples)
        tmp2 = np.ones(n_samples) * C
        h = cvxopt.matrix(np.hstack((tmp1, tmp2)))

    solution = cvxopt.solvers.qp(P, q, G, h, A, b)

    alphas = np.ravel(solution['x'])

수식에 대한 설명만 보고 바로 코드로 넘어오면 혼란스럽기 그지 없습니다. cvxopt 라이브러리를 사용해야하기 때문에 많은 부분에 변화가 생겼기 때문입니다.

우선 cvxopt.solver.oq 가 어떤 식으로 최적화를 하는지 보겠습니다. cvxopt.solvers.qp(P, q, G, h, A, b) 가 직접적으로 $\alpha$ 를 도출하는 명령입니다. 이 스크립트의 연산과정을 인간의 눈으로 볼 수 있게 바꾸면 다음과 같습니다.

$\min \frac {1} {2} x^{T}Px -q^{T}x$ $s.t. \; Gx \ge h \; \; and \;\; Ax=b$

다른 방법이 없기 때문에 울며 겨자먹기로 위의 형태에 알맞게 앞서 구했던 Dual problem 을 아래와 같이 고쳐줍니다.

$\min \frac {1} {2} \alpha^{T}K \alpha -1^{T}\alpha$ $s.t. \; \alpha_i \le 0 \;\; and \;\; y^{T}\alpha=0$

코드를 자세히 살펴보면 C에 따라서 G와 h를 구하는 방법이 다르게 나와있습니다. 이 부분은 바로 뒤에 나올 C-SVM 에 관한 내용입니다.

cvxopt 를 통해 도출해낸 alpha 에서 1e-5 이하인 값은 걸러냅니다.
이보다 더 작은 값은 너무 작아 0이나 마찬가지 이므로
차후 w나 b를 도출하는 데에 있어 노이즈가 될 수 있기 때문입니다.
아래에 해당 내용이 있습니다.

 # Support vectors have non zero alphas
    cdt = alphas > 1e-5
    idx = np.where(cdt)[0]


    self.alphas = alphas[cdt]
    self.sv_x = X[cdt]
    self.sv_y = y[cdt]
    print("%d support vectors among %d points" % (len(self.alphas), n_samples))

이제 다음로 넘어가서 w 와 b를 구하는 과정에 대해서 보도록 하겠습니다.

  # b = y - summation(alphas*y*k)
    self.b = 0
    for n in range(len(self.alphas)):
        self.b += self.sv_y[n] - np.sum(self.alphas * self.sv_y * K[idx[n],cdt])
    self.b /= len(self.alphas)

    # Weight vector, w = sum(alph*y*x)
    if self.kernel == linear_kernel:
        self.w = np.zeros(n_features)
        for n in range(len(self.alphas)):
            self.w += self.alphas[n] * self.sv_y[n] * self.sv_x[n]

    else:
        self.w = None

C-SVM

앞서 봤던 SVM은 마진 사이에 관측치가 들어올 수 없는 strict한 hard margin 이었습니다.
하지만 어느 정도 ‘감당할 수 있는 x(관측치)’가 마진 사이에 있도록 허락하면
더 큰 Margin을 얻을 수 있지 않을까 라는 아이디어를 통해 C-SVM 이라는 알고리즘이 탄생했습니다.
아래의 그림을 통해 C-SVM이 어떻게 작동하는지 감을 잡을 수 있을 것입니다.

C-svm

위의 그림에서 점선을 벗어난 파란색 점과 빨간색 점을 각각 하나씩 볼 수 있습니다. 저렇게 해당 관측치가 margin 을 벗어난 정도를 $\xi$ 로 나타내게 되고, 목적함수에 그대로 반영됩니다.

C-SVM 의 목적함수는 다음과 같습니다.

$\min { \frac { 1 }{ 2 } { \left\| w \right\| }_{ 2 }^{ 2 } } +C\sum _{ i=1 }^{ n }{ { \xi }_{ i } }$

위의 목적식을 두 부분으로 나눠서 살펴보겠습니다.

$\min { \frac { 1 }{ 2 } { \left| w \right| }_{ 2 }^{ 2 } }$ 은 마진의 역수에 대한 최소화이므로, 마진을 최대화하는 것과 같습니다.

$\min C\sum _{i=1}^{n}{\xi_i}$ 은 C-svm 이 Margin을 벗어나는 것을 허용하지만 그로 인해 생기는 페널티를 최소화하고자 하는 것을 알 수 있게 해주는 부분입니다.

$s.t. \quad { y }_{ i }({ w }^{ T }{ x }_{ i }+b)\ge 1-{ \xi }_{ i },\quad { \xi }_{ i }\ge 0$

위의 제약식 역시 새로 생겨난 penalty term 인 $\xi_i$ 에의한 변화 외에는 기존의 svm과 다른 점이 없습니다.
따라서 이 식에 대한 Dual problem 화는 생략하고 바로 결과물을 가져오도록 하겠습니다.

아래의 식이 C-svm 의 Dual 식입니다.

$L_d$ 의 제약식은 다음과 같습니다.

$\sum _{ i=1 }^{ n }{ { \alpha }_{ i }{ y }_{ i } } =0 \\0\le { \alpha }_{ i }\le C,\quad i=1,...,n$

C-SVM 결과

KKT condition 에 의해 다음을 알 수 있습니다.

$\alpha \left( y_i \left(w^{T}x_i + b \right)-1+\xi_{i} \right) = 0$

위 식에 따라 support vector 는 $\alpha \ne 0$ 일 때임을 알 수 있습니다.

또한 주어진 아래의 식을 통해서 C-SVM를 세 가지 케이스로 나눌 수 있게 됩니다.

$C- \alpha_{i}-u_{i} = 0 \; and \; u_{i}\xi_i = 0$

case 1. $\alpha_i = 0$ 이면 관측치는 non-support vector

$\alpha_i = 0$ 이면 $y_i \left( w^{T}x_i +b \right) > 1-\xi_i$ 이어야 합니다.
그러면 $u_i$ 는 $C$ 와 같아야 하고, $x_i=0$이 됩니다.
결국 $y_i \left( w^Tx+b \right) > 1$ 이라는 결과가 나오므로 해당 관측치는 support vector가 될 수 없습니다.

case 2. $0 < \alpha_i < C$ 면 관측치는 margin 위의 Suppor Vector

이 경우 $u_i > 0$ 이므로, $\xi_i$ 는 0이 됩니다.
따라서 $y_i \left( w^Tx+b \right) = 1$ 이므로, 해당 데이터는 support vector 가 됩니다.

case 3. $\alpha_i = C$ 이면 관측치는 마진 밖의 vector

당 케이스의 경우 $u_i = 0$ 이므로, $\xi_i > 0$ 이 됩니다.
$\xi$ 는 마진을 벗어난 정도를 나타내므로 0보다 크면 해당 데이터는 마진을 벗어난 값입니다.

C의 역할

C 는 Regularization을 해줍니다.
Large C 의 경우에는 좁은 마진과 적은 수의 support vector 가 생성됩니다.
Small C 의 경우에는 넓은 마진과 많은 수의 support vector 가 생성됩니다.

그 이유는 C가 클 수록 목적식에서 $\xi_i$ 가 미치는 영향이 커지게 되고, 그에 따라 마진이 적어지게 됩니다.
반대로 C가 작으면 $\xi_i$ 가 적당히 크더라도 그 영향이 작아지므로 마진에 미치는 영향이 작아질 수 밖에 없습니다.
C에 의한 영향은 아래의 그림을 보면 명확히 다가올 것입니다.

Kernel-SVM

선형성과 내적으로 분류문제를 해결하는 SVM 이 해결하지 못하는 문제가 있습니다. 아래의 그래프와 같은 XOR problem은 하나의 직선으로 절대 분류할 수 없는 문제입니다. Kernel SVM 은 이렇게 선형으로 풀 수 없는 문제를 Kernel 이라는 트릭을 통해서 해결하려 합니다.

‘본래의 차원에서는 도저히 파란색과 빨간색 사이를 가를 수 없지만, 똑같은 데이터를 더 높은 차원으로 확장시킨다면 자를 수 있지 않을까?’라는 것이 Kernel SVM 의 메인 아이디어입니다. 아래의 그림이 그것을 보여주는 정확한 예입니다.

출처:https://quantdare.com/svm-versus-a-monkey/

출처:https://quantdare.com/svm-versus-a-monkey/
이러한 문제를 해결하기 위한 목적식과 제약식은 C-SVM의 것과 크게 다르지 않습니다. $$ \min \frac{1}{2} \parallel w \parallel + C \sum_{i=1}^N \xi_i $$ $$ s.t. \quad y_i \left( w^T \Phi \left(x_i \right) + b \right) \ge 1-\xi_i, \;\; \xi_i \ge 0, \forall i $$ 이전의 것과 다른 점은 제약식에 $x_i$ 를 감싸는 function $\Phi$가 생겼다는 것입니다. 이 $\Phi$가 저차원의 데이터를 고차원으로 변형해주는 mapping function 입니다.
바로 Dual 문제로 넘어가겠습니다. $$ max L_D = \sum_{i=1}^n \alpha_i - \frac{1}{2} \sum_{i=1}^n \sum_{j=1}^n {\alpha_i \alpha_j y_i y_j \Phi \left(x_i \right)^T \Phi \left(x_j \right) } $$ $$ s.t. \quad \sum_{i=1}^n =0 \; and \; 0 \le \alpha_i \le C $$ python code ----- SVM 실행 위한 전체코드는 다음과 같습니다. 해당 코드는 어떤 깃허브의 코드를 토대로 약간의 변경이 있습니다만 출처를 기억하지 못해 첨부하지 못했습니다. 속히 찾아서 반영토록하겠습니다. $$ a $$ ```python import numpy as np from numpy import linalg import cvxopt import cvxopt.solvers import pylab as pl ``` ```python def linear_kernel(x1, x2): return np.dot(x1, x2) def poly_kernel(x, y, p=4): return (np.dot(x, y)+1) ** p def RBF_kernel(x, y, sigma=5.0): return np.exp(-linalg.norm(x-y)**2 / (2 * (sigma ** 2))) class SVM: def __init__(self, kernel=RBF_kernel, C=None): self.kernel = kernel self.C = C if self.C is not None: self.C = float(self.C) def fit(self, X, y): n_samples, n_features = X.shape # Gram matrix K = np.zeros((n_samples, n_samples)) for i in range(n_samples): for j in range(n_samples): K[i,j] = self.kernel(X[i], X[j]) P = cvxopt.matrix(np.outer(y,y) * K) q = cvxopt.matrix(np.ones(n_samples) * -1) A = cvxopt.matrix(y, (1,n_samples)) b = cvxopt.matrix(np.zeros(1)) if self.C is None: G = cvxopt.matrix(np.diag(np.ones(n_samples) * -1)) h = cvxopt.matrix(np.zeros(n_samples)) else: tmp1 = np.diag(np.ones(n_samples) * -1) tmp2 = np.identity(n_samples) G = cvxopt.matrix(np.vstack((tmp1, tmp2))) tmp1 = np.zeros(n_samples) tmp2 = np.ones(n_samples) * C h = cvxopt.matrix(np.hstack((tmp1, tmp2))) solution = cvxopt.solvers.qp(P, q, G, h, A, b) alphas = np.ravel(solution['x']) # Support vectors have non zero alphas cdt = alphas > 1e-5 idx = np.where(cdt)[0] self.alphas = alphas[cdt] self.sv_x = X[cdt] self.sv_y = y[cdt] print("%d support vectors among %d points" % (len(self.alphas), n_samples)) # b = y - summation(alphas*y*k) self.b = 0 for n in range(len(self.alphas)): self.b += self.sv_y[n] - np.sum(self.alphas * self.sv_y * K[idx[n],cdt]) self.b /= len(self.alphas) # Weight vector, w = sum(alph*y*x) if self.kernel == linear_kernel: self.w = np.zeros(n_features) for n in range(len(self.alphas)): self.w += self.alphas[n] * self.sv_y[n] * self.sv_x[n] else: self.w = None ## project => separate def separate(self, X): if self.w is not None: return np.dot(X, self.w) + self.b else: y_predict = np.zeros(len(X)) for i in range(len(X)): s = 0 for a, y, x in zip(self.alphas, self.sv_y, self.sv_x): s += a * y * self.kernel(X[i], x) y_predict[i] = s return y_predict + self.b def predict(self, X): return np.sign(self.separate(X)) def plot_margin(X1_train, X2_train, svm): def f(x, w, b, c=0): # given x, return y such that [x,y] in on the line # w.x + b = c return (-w[0] * x - b + c) / w[1] pl.plot(X1_train[:,0], X1_train[:,1], "ro") pl.plot(X2_train[:,0], X2_train[:,1], "bo") # support vector pl.scatter(svm.sv_x[:,0], svm.sv_x[:,1], s=100, c="g") # w.x + b = 0 , 중간선 a0 = -5 # 임의의 수를 넣음 a1 = f(a0, svm.w,svm.b) # (-w[0] * x - b + c) / w[1] b0 = 5 # 임의의 수를 넣음 b1 = f(b0, svm.w, svm.b) # (-w[0] * x - b + c) / w[1] # slope = (a1-b1)/(a0-b0) # (a0,a1) => (b1,b2) pl.plot([a0,b0], [a1,b1], "k") # w.x + b = 1 a0 = -4 a1 = f(a0, svm.w, svm.b, 1) b0 = 4 b1 = f(b0, svm.w, svm.b, 1) pl.plot([a0,b0], [a1,b1], "k--") # w.x + b = -1 a0 = -4; a1 = f(a0, svm.w, svm.b, -1) b0 = 4; b1 = f(b0, svm.w, svm.b, -1) pl.plot([a0,b0], [a1,b1], "k--") pl.axis("tight") pl.show() # clf => SVM def plot_contour(X1_train, X2_train, model): pl.plot(X1_train[:, 0], X1_train[:, 1], "ro") pl.plot(X2_train[:, 0], X2_train[:, 1], "bo") pl.scatter(model.sv_x[:, 0], model.sv_x[:, 1], s=100, c="g") X1, X2 = np.meshgrid(np.linspace(-6, 6, 50), np.linspace(-6, 6, 50)) X = np.array([[x1, x2] for x1, x2 in zip(np.ravel(X1), np.ravel(X2))]) Z = model.separate(X).reshape(X1.shape) pl.contour(X1, X2, Z, [0.0], colors='k', linewidths=1, origin='lower') pl.contour(X1, X2, Z + 1, [0.0], colors='grey', linewidths=1, origin='lower') pl.contour(X1, X2, Z - 1, [0.0], colors='grey', linewidths=1, origin='lower') pl.axis("tight") pl.show() def test_linear(): X1, y1, X2, y2 = gen_lin_separable_data() X_train, y_train = split_train(X1, y1, X2, y2) X_test, y_test = split_test(X1, y1, X2, y2) clf = SVM(kernel=linear_kernel) clf.fit(X_train, y_train) y_predict = clf.predict(X_test) correct = np.sum(y_predict == y_test) print("%d out of %d predictions correct" % (correct, len(y_predict))) plot_margin(X_train[y_train==1], X_train[y_train==-1], clf) def test_non_linear(): X1, y1, X2, y2 = gen_non_lin_separable_data() X_train, y_train = split_train(X1, y1, X2, y2) X_test, y_test = split_test(X1, y1, X2, y2) svm = SVM(gaussian_kernel) svm.fit(X_train, y_train) y_predict = clf.predict(X_test) correct = np.sum(y_predict == y_test) print("%d out of %d predictions correct" % (correct, len(y_predict))) plot_contour(X_train[y_train==1], X_train[y_train==-1], svm) def test_soft(): X1, y1, X2, y2 = gen_lin_separable_overlap_data() X_train, y_train = split_train(X1, y1, X2, y2) X_test, y_test = split_test(X1, y1, X2, y2) svm = SVM(kernel=linear_kernel,C=0.1) svm.fit(X_train, y_train) y_predict = svm.predict(X_test) correct = np.sum(y_predict == y_test) print("%d out of %d predictions correct" % (correct, len(y_predict))) plot_margin(X_train[y_train==1], X_train[y_train==-1], svm) def gen_non_lin_separable_data(): mean1 = [-1, 2] mean2 = [1, -1] mean3 = [4, -4] mean4 = [-4, 4] cov = [[1.0,0.8], [0.8, 1.0]] X1 = np.random.multivariate_normal(mean1, cov, 50) X1 = np.vstack((X1, np.random.multivariate_normal(mean3, cov, 50))) y1 = np.ones(len(X1)) X2 = np.random.multivariate_normal(mean2, cov, 50) X2 = np.vstack((X2, np.random.multivariate_normal(mean4, cov, 50))) y2 = np.ones(len(X2)) * -1 return X1, y1, X2, y2 def split_train(X1, y1, X2, y2): X1_train = X1[:90] y1_train = y1[:90] X2_train = X2[:90] y2_train = y2[:90] X_train = np.vstack((X1_train, X2_train)) y_train = np.hstack((y1_train, y2_train)) return X_train, y_train def split_test(X1, y1, X2, y2): X1_test = X1[90:] y1_test = y1[90:] X2_test = X2[90:] y2_test = y2[90:] X_test = np.vstack((X1_test, X2_test)) y_test = np.hstack((y1_test, y2_test)) return X_test, y_test ```

Jeong Minsung