機器學習 支持向量機（SVM） 從理論到放棄，從代碼到理解
基本概念
支持向量機（support vector machines，SVM）是一種二分類模型，它的基本模型是定義在特徵空間上的間隔最大的線性分類器。支持向量機還包括核技巧，這使它成爲實質上的非線性分類器。其學習策略就是間隔最大化，可形式化爲一個求解凸二次規劃（convex quadratic programming）的問題，也等價於正則化的合頁損失函數的最小化問題，支持向量機的學習算法是求解凸二次規劃的最優化算法。html
支持向量機學習方法包含構建由簡至繁的模型：線性可分支持向量機（linear support vector machine in linearly separable case）、線性支持向量機（linear support vector machine）及非線性支持向量機（non-linear support vector machine）。當訓練數據線性可分時，經過硬間隔最大化（hard margin maximization），學習一個線性的分類器；當訓練數據近似線性可分時，經過軟間隔最大化（soft margin maximization），也學習一個線性的分類器；當訓練數據不可分時，經過核技巧（kernel trick）及軟間隔最大化，訓練非線性的分類器。python
最大邊緣超平面
對於一個二分類問題，若數據集是線性可分的，那麼咱們能夠找到這樣一個超平面，使得數據的兩個label分別位於平面兩側。而且，能夠注意到，這樣的超平面咱們能夠找到無數個。算法
然而，雖然有無數個這樣的超平面，咱們並不能保證這些超平面在未知實例上的效果都同樣的好。所以，咱們須要找到一個具備很好的泛化偏差的超平面。這裏，咱們能夠選取最大邊緣超平面。下面將講述緣由。數組
函數間隔
下圖中有A，B，C三點，表示3個實例，對於這3個實例的預測，咱們預測的確信度會是A>B>C這樣一個結果。通常來講，一個點距離超平面的遠近能夠表示分類預測的確信程度。在超平面 $w \cdot x+b=0$ 肯定的狀況下， $|w \cdot x +b|$ 可以相對地表示點 $x$ 距離超平面地遠近。而 $w \cdot x + b$ 的符號與相似標記 $y$ 的符號是否一致可以表示分類是否正確。因此可用量 $y(w \cdot x + b)$ 來表示分類的正確性及確信度，這就是函數間隔（function margin）。緩存
幾何間隔
函數間隔能夠表示分類預測的正確性及確信度，可是選擇分離超平面時，只有函數間隔還不夠，由於只要成比例地該邊$w和b$，超平面並無改變，但函數間隔卻變爲原來地兩倍。所以，咱們須要對分離超平面的法向量 $w$ 加某些約束，如規範化。這時，函數間隔便成了幾何間隔（geometric margin）。app
間隔最大化
對訓練數據集找到幾何間隔最大的超平面意味着以充分大的確信度對訓練數據進行分類。即，不只將正負實例點分開，並且對最難分的實例點也有足夠大的確信度將它們分開，這樣的超平面對未知的新實例有很好的分類預測能力。這也是爲何咱們要選取最大邊緣超平面的緣由。dom
線性可分支持向量機
最大間隔法
最優化
運用拉格朗日對偶性，經過求解對偶問題（dual problem）獲得原始問題的最優解。機器學習
優勢：ide
- 對偶問題每每更易求解
- 天然引入核函數，進而推廣到非線性分類問題
首先構建拉格朗日函數，對每一個不等式約束引進拉格朗日乘子 $\alpha_i \ge 0,i=1,2, \cdots,N$ ,定義拉格朗日函數： $$ L(w,b,\alpha)=\frac{1}{2}||w||^2-\sum_{i=1}^{N}\alpha_iy_i(w\cdot x_i+b)+\sum_{i=1}^{N}\alpha_i $$ 其中，$\alpha =(\alpha_1,\alpha_2, \cdots ,\alpha_N)^T$ 爲拉格朗日乘子向量。函數
根據拉格朗日對偶性，原始問題的對偶問題就是極大極小問題： $$ \max_a\min_{w,b}L(w,b,a) $$ 附上推導過程
線性支持向量機
當數據集沒法線性可分時，以下圖所示
這時，咱們能夠對每個樣本點引進一個鬆弛變量 $\zeta_i \ge 0 $ ，使函數間隔加上鬆弛變量大於等於1，這樣，約束條件就變成 $$ y_i(w \cdot x+b)\ge 1-\zeta_i $$ 同時，對每一個鬆弛變量 $ \zeta_i $ ，支付一個代價 $\zeta_i $ 。目標函數由原來的$\frac{1}{2}||w||^2$ 變成 $$ \frac{1}{2}||w||^2+C\sum_{i=1}^{N}\zeta_i $$ 這裏，$ C > 0 $ 稱爲懲罰參數，通常由應用問題決定。最小化目標函數包含兩層含義：使 $ \frac{1}{2}||w||^2 $ 儘可能小，即間隔儘可能大，同時使偏差分類點的個數儘可能小，C是調和兩者的係數。
同上述推導方法，咱們可寫出下面算法
非線性支持向量機
有時，咱們面對的是非線性問題，以下圖左。非線性問題每每很差求解，因此但願能用解線性分類問題的方法解決這個問題，採起的方法是進行一個非線性變換，將非線性問題變換爲線性問題，經過解變換後的線性問題的方法求解原來的非線性問題。
原空間： $$ \mathcal{X} \subset \mathbf{R}^{2}, x=\left(x^{(1)}, x^{(2)}\right)^{\mathrm{T}} \in \mathcal{X} $$ 新空間： $$ \mathcal{Z} \subset \mathbf{R}^{2}, z=\left(z^{(1)}, z^{(2)}\right)^{\mathrm{T}} \in \mathcal{Z} \quad z=\phi(x)=\left(\left(x^{(1)}\right)^{2},\left(x^{(2)}\right)^{2}\right)^{\mathrm{T}} $$
核函數
設 $\mathcal{X}$ 是輸入空間(歐氏空間 $\mathbf{R}^{n}$ 的子集或離散集合)，又設 $\mathcal{H}$ 爲特徵空間(希爾伯特空間)，若是存在一個從 $\mathcal{X}$ 到 $\mathcal{H}$ 的映射 $$ \phi(x) : \mathcal{X} \rightarrow \mathcal{H} $$ 使得對全部 $$ x, z \in \mathcal{X} $$ 函數K(x,z)知足條件 $$ K(x, z)=\phi(x) \cdot \phi(z) $$ 則稱 K(x,z) 爲核函數， $\phi(x) $ 爲映射函數， 式中 $\phi(x) \cdot \phi(z)$ 爲 $\phi(x) $ 和 $\phi(z)$ 的內積。
在學習與預測中只定義核函數K(x,z)，而不顯式地定義映射函數，一般，直接計算K(x,z)比較容易，而經過 $\phi(x) $ 和 $\phi(z)$ 計算K(x, z)並不容易。
注：φ是輸入空間 $\mathbf{R}^{n}$ 到特徵空間 $\mathcal{H}$ 的映射，特徵空間 $\mathcal{H}$ 通常是高維，映射能夠不一樣。
引入核函數後，目標函數改成 $$ W(\alpha)=\frac{1}{2} \sum_{i=1}^{N} \sum_{j=1}^{N} \alpha_{i} \alpha_{j} y_{i} y_{j} K\left(x_{i}, x_{j}\right)-\sum_{i=1}^{N} \alpha_{i} $$
正定核
水平有限，不會推導，就直接給結論好了。。orz
一般咱們所說的核函數就是正定核函數，設 $\mathrm{K} : \mathcal{X} \times \mathcal{X} \rightarrow \mathbf{R}$ 是對稱函數，則K(x,z)爲正定核函數 的充要條件是對任意 $x_{i} \in \mathcal{X}, \quad i=1,2, \cdots, m$ ，K(x,z)對應的Gram矩陣 $$ K=\left[K\left(x_{i}, x_{j}\right)\right]_{m \times m} $$ 是半正定的。
這必定義在構造核函數時頗有用。但對於一個具體函數 K(x,z) 來講，檢驗它是否爲正定核函數並不容易，由於 要求對任意有限輸入集 $\left{x_{1}, x_{2}, \cdots, x_{m}\right}$ 驗證K對應的 Gram矩陣是否爲半正定的。在實際問題中每每應用己有的核函數。
經常使用核函數
- 多項式核函數（Polynomial kernel function)
$$ K(x, z)=(x \cdot z+1)^{p} $$
對應的支持向量機爲P次多項式分類器，分類決策函數： $$ f(x)=\operatorname{sign}\left(\sum_{i=1}^{N_{s}} a_{i}^{} y_{i}\left(x_{i} \cdot x+1\right)^{p}+b^{}\right) $$
- 高斯核函數 （Gaussian Kernel Function)
$$ K(x, z)=\exp \left(-\frac{|x-z|^{2}}{2 \sigma^{2}}\right) $$
決策函數： $$ f(x)=\operatorname{sign}\left(\sum_{i=1}^{N_{s}} a_{i}^{} y_{i} \exp \left(-\frac{|x-z|^{2}}{2 \sigma^{2}}\right)+b^{}\right) $$
- 字符串核函數（string kernel function）：
核函數不只能夠定義在歐式空間上，還能夠定義在離散數據的集合上。好比，字符串核是定義在字符串集合上的核函數。字符串核函數在文本分類、信息檢索、生物信息學等方面都有應用。 $$ k_n(s,t)=\sum_{u\in \sum^n}[\phi_n(s)]u[\phi_n(t)]u=\sum{u\in \sum^n}\sum{(i,j):s(i)=t(j)=u}\lambda^{l(i)}\lambda^{l(j)} $$
字符串核函數 $k_n(s,t)$ 給出了字符串s和t中長度等於n的全部字串組成的特徵向量的餘弦類似度。
算法
序列最小優化（SMO）算法
基本思路
若是全部變量的解都知足此最優化問題的KKT條件，那麼獲得解； 不然，選擇兩個變量，固定其它變量，針對這兩個變量構建一個二次規劃問題，稱爲子問題，可經過解析方法求解，提升計算速度。
子問題的兩個變量：一個是違反KKT條件最嚴重的那個，另 一個由約束條件自動肯定。假設 $\alpha_1,\alpha_2$ 爲兩個變量， $\alpha_3,\alpha_4, \cdots , \alpha_ N$ 固定，那麼由等式約束可知 $$ \alpha_1=-y_1\sum_{i=2}^{N}\alpha_iy_i $$
變量選擇
SMO算法在每一個子問題中選擇兩個變量優化，其中至少一個變量是違反KKT條件的。
- 第一個變量的選擇：外循環，違反KKT最嚴重的樣本點，檢驗樣本點是否知足KKT條件：
$$ \alpha_{i}=0 \Leftrightarrow y_{i} g\left(x_{i}\right) \geqslant 1 \ \begin{aligned} 0<\alpha_{i} &<C \Leftrightarrow y_{i} g\left(x_{i}\right)=1 \ \alpha_{i} &=C \Leftrightarrow y_{i} g\left(x_{i}\right) \leqslant 1 \ g\left(x_{i}\right) &=\sum_{j=1}^{N} \alpha_{j} y_{j} K\left(x_{i}, x_{j}\right)+b \end{aligned} $$
-
第二個變量的檢查: 內循環， 選擇的標準是但願能使目標函數有足夠大的變化，即對應 $|E_1-E_2|$ 最大 ,若是內循環經過上述方法找到的點不能使目標函數有足夠的降低，則：遍歷間隔邊界上的樣本點，測試目標函數降低 ，若是降低不大，則遍歷全部樣本點 ，若是依然降低不大，則丟棄外循環點，從新選擇。
-
每次完成兩個變量的優化後，從新計算 $b,E_i$
算法
代碼
這是我用的數據集數據的分佈
感受仍是看代碼理解的快，先貼一份根據《機器學習實戰》上寫的代碼
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import random
import math
from sklearn.model_selection import train_test_split
from sklearn.svm import SVC
"""
Author:
ITryagain
Modify:
2019-04-09
"""
def loadDataSet(fileName):
"""
讀取數據
Parameters:
fileName - 文件名
Returns:
dataMat - 數據矩陣
labelMat - 數據標籤
"""
dataMat = []
labelMat = []
fr = open(fileName)
for line in fr.readlines(): # 逐行讀取，濾除空格等
lineArr = line.strip().split('\t')
dataMat.append([float(lineArr[0]), float(lineArr[1])]) # 添加數據
labelMat.append(float(lineArr[2])) # 添加標籤
return np.array(dataMat), np.array(labelMat)
class optStruct:
def __init__(self, dataMatIn, classLabels, C, toler, kTup):
"""
:param dataMatIn: 輸入數據 [X1, X2, ... , XN]
:param classLabels: 分類標籤 [y]
:param C: 鬆弛變量
:param toler: 容錯率
"""
self.X = dataMatIn
self.labelMat = classLabels
self.C = C
self.tol = toler
self.m = np.shape(dataMatIn)[0]
self.alphas = np.mat(np.zeros((self.m, 1)))
self.b = 0
self.eCache = np.mat(np.zeros((self.m, 2))) # 偏差緩存
self.k = np.mat(np.zeros((self.m, self.m)))
for i in range(self.m):
self.k[:, i] = kernelTrans(self.X, self.X[i, :], kTup)
def calcEk(oS, k):
"""
:param oS:
:param k: 第k行
:return: 預測值與實際值的差
"""
fXk = float(np.multiply(oS.alphas, oS.labelMat).T * oS.k[:, k]) + oS.b
Ek = fXk - float(oS.labelMat[k])
return Ek
def clipAlpha(aj, H, L):
"""
調整aj的值，使aj處於 L<=aj<=H
:param aj: 目標值
:param H: 最大值
:param L: 最小值
:return:
aj 目標值
"""
if aj > H:
aj = H
elif L > aj:
aj = L
return aj
def selectJrand(i, m):
"""
隨機選擇一個整數
:param i: 第一個alpha的下標
:param m: 全部alpha的數目
:return:
j 返回一個不爲i的隨機數，在0~m之間的整數值
"""
j = i
while j == i:
j = int(random.uniform(0, m))
return j
def selectJ(i, oS, Ei):
"""
內循環的啓發式方法，選擇第二個alpha的值
:param i: 第一個alpha的下標
:param oS:
:param Ei: 預測結果與真實結果比對，計算偏差Ei
:return:
j 隨機選出的第j一行
Ej 預測結果與真實結果比對，計算偏差Ej
"""
maxK = -1
maxDeltaE = 0
Ej = 0
oS.eCache[i] = [1, Ei] # 將輸入值Ei在緩存中設置成爲有效的。這裏的有效意味着它已經計算好了。
validEcacheList = np.nonzero(oS.eCache[:, 0].A)[0] # 非零E值的行的list列表，所對應的alpha值
if(len(validEcacheList)) > 1:
for k in validEcacheList: # 在全部的值上進行循環，並選擇其中使得改變最大的那個值
if k ==1:
continue
Ek = calcEk(oS, k)
deltaE = abs(Ei - Ek)
if(deltaE > maxDeltaE):
maxk = k
maxDeltaE = deltaE
Ej = Ek
return maxK, Ej
else: # 若是是第一次循環，則隨機選擇一個alpha值
j = selectJrand(i, oS.m)
Ej = calcEk(oS, j)
return j, Ej
def updateEk(oS, k):
"""
計算偏差值並存入緩存中
:param oS:
:param k: 某一列的行號
"""
Ek = calcEk(oS, k)
oS.eCache[k] = [1, Ek]
def kernelTrans(X, A, kTup):
"""
核轉換函數
:param X: dataMatIn數據集
:param A: dataMatIn數據集的第i行的數據
:param kTup:核函數的信息
:return:
"""
m, n = np.shape(X)
K = np.mat(np.zeros((m, 1)))
if kTup[0] == 'lin':
# linear kernel: m*n * n*1 = m*1
K = X * A.T
elif kTup[0] == 'rbf':
for j in range(m):
deltaRow = X[j, :] - A
K[j] = deltaRow * deltaRow.T
# 徑向基函數的高斯版本
# print(len(K))
# K = math.exp(K / (-1 * kTup[1] ** 2))
for i in range(m):
K[i] = math.exp(K[i] / (-1 * kTup[1] ** 2)) # divide in NumPy is element-wise not matrix like Matlab
else:
raise NameError('Houston We Have a Problem -- That Kernel is not recognized')
return K
def innerL(i, oS):
"""
內循環代碼
:param i:具體的某一行
:param oS:
:returns:
0 找不到最優的值
1 找到了最優的值，而且oS.Cache到緩存中
"""
Ei = calcEk(oS, i)
# 約束條件 (KKT條件是解決最優化問題的時用到的一種方法。咱們這裏提到的最優化問題一般是指對於給定的某一函數，求其在指定做用域上的全局最小值)
# 0<=alphas[i]<=C，但因爲0和C是邊界值，咱們沒法進行優化，由於須要增長一個alphas和下降一個alphas。
# 表示發生錯誤的機率：labelMat[i]*Ei 若是超出了 toler， 才須要優化。至於正負號，咱們考慮絕對值就對了。
'''
# 檢驗訓練樣本(xi, yi)是否知足KKT條件
yi*f(i) >= 1 and alpha = 0 (outside the boundary)
yi*f(i) == 1 and 0<alpha< C (on the boundary)
yi*f(i) <= 1 and alpha = C (between the boundary)
'''
if((oS.labelMat[i] * Ei < -oS.tol) and (oS.alphas[i] < oS.C)) or ((oS.labelMat[i] * Ei > oS.tol) and (oS.alphas[i] > 0)):
# 選擇最大的偏差對應的j進行優化。效果更明顯
j ,Ej = selectJ(i, oS, Ei)
alphaIold = oS.alphas[i].copy()
alphaJold = oS.alphas[j].copy()
# L和H用於將alphas[j]調整到0-C之間。若是L==H，就不作任何改變，直接return 0
if(oS.labelMat[i] != oS.labelMat[j]):
L = max(0, oS.alphas[j] - oS.alphas[i])
H = min(oS.C, oS.C + oS.alphas[j] - oS.alphas[i])
else:
L = max(0, oS.alphas[j] + oS.alphas[i] - oS.C)
H = min(oS.C, oS.alphas[j] + oS.alphas[i])
if L ==H:
print("L==H")
return 0
# eta是alphas[j]的最優修改量，若是eta==0，須要退出for循環的當前迭代過程
# 參考《統計學習方法》李航-P125~P128<序列最小最優化算法>
eta = 2.0 * oS.k[i,j] - oS.k[i, i] - oS.k[j, j]
if eta >= 0:
print("eta>=0")
return 0
# 計算出一個新的alphas[j]值
oS.alphas[j] -= oS.labelMat[j] * (Ei - Ej) / eta
# 並使用輔助函數，以及L和H對其進行調整
oS.alphas[j] = clipAlpha(oS.alphas[j], H, L)
# 更新偏差緩存
updateEk(oS, j)
# 檢查alpha[j]是否只是輕微的改變，若是是的話，就退出for循環。
if abs(oS.alphas[j] - alphaJold) < 0.00001:
print("j not moving enough")
return 0
# 而後alphas[i]和alphas[j]一樣進行改變，雖然改變的大小同樣，可是改變的方向正好相反
oS.alphas[i] += oS.labelMat[j] * oS.labelMat[i] * (alphaJold - oS.alphas[j])
# 更新偏差緩存
updateEk(oS, i)
# 在對alpha[i], alpha[j] 進行優化以後，給這兩個alpha值設置一個常數b。
# w= Σ[1~n] ai*yi*xi => b = yj Σ[1~n] ai*yi(xi*xj)
# 因此： b1 - b = (y1-y) - Σ[1~n] yi*(a1-a)*(xi*x1)
# 爲何減2遍？ 由於是 減去Σ[1~n]，正好2個變量i和j，因此減2遍
b1 = oS.b - Ei - oS.labelMat[i] * (oS.alphas[i] - alphaIold) * oS.k[i, i] - oS.labelMat[j] * (
oS.alphas[j] - alphaJold) * oS.k[i, j]
b2 = oS.b - Ej - oS.labelMat[i] * (oS.alphas[i] - alphaIold) * oS.k[i, j] - oS.labelMat[j] * (
oS.alphas[j] - alphaJold) * oS.k[j, j]
if (0 < oS.alphas[i]) and (oS.C > oS.alphas[i]):
oS.b = b1
elif (0 < oS.alphas[j]) and (oS.C > oS.alphas[j]):
oS.b = b2
else:
oS.b = (b1 + b2) / 2.0
return 1
else:
return 0
def smoP(dataMatIn, classLabels, C, toler, maxIter, kTup = ('lin', 0)):
"""
SMO算法外循環
:param dataMatIn: 輸入數據集
:param classLabels: 標籤
:param C: 鬆弛變量(常量值)，容許有些數據點能夠處於分隔面的錯誤一側。
控制最大化間隔和保證大部分的函數間隔小於1.0這兩個目標的權重。
能夠經過調節該參數達到不一樣的結果。
:param toler: 容錯率
:param maxIter: 退出前最大的循環次數
:param kTup: 核函數
:return:
"""
# 建立一個 optStruct 對象
oS = optStruct(np.mat(dataMatIn), np.mat(classLabels).transpose(), C, toler, kTup)
iter = 0
entireSet = True
alphaPairsChanged = 0
# 循環遍歷：循環maxIter次 而且 （alphaPairsChanged存在能夠改變 or 全部行遍歷一遍）
# 循環迭代結束 或者 循環遍歷全部alpha後，alphaPairs仍是沒變化
while iter < maxIter and (alphaPairsChanged > 0 or entireSet):
alphaPairsChanged = 0
# 當entireSet=true or 非邊界alpha對沒有了；就開始尋找 alpha對，而後決定是否要進行else。
if entireSet:
# 在數據集上遍歷全部可能的alpha
for i in range(oS.m):
# 是否存在alpha對，存在就+1
alphaPairsChanged += innerL(i, oS)
print("fullSet, iter: %d i:%d, pairs changed %d" % (iter, i, alphaPairsChanged))
iter += 1
# 對已存在 alpha對，選出非邊界的alpha值，進行優化。
else:
# 遍歷全部的非邊界alpha值，也就是不在邊界0或C上的值。
nonBoundIs = np.nonzero((oS.alphas.A > 0) * (oS.alphas.A < C))[0]
for i in nonBoundIs:
alphaPairsChanged += innerL(i, oS)
print("non_bound, iter: %d i:%d, pairs changed %d" % (iter, i, alphaPairsChanged))
iter += 1
# 若是找到alpha對，就優化非邊界alpha值，不然，就從新進行尋找，若是尋找一遍 遍歷全部的行仍是沒找到，就退出循環。
if entireSet:
entireSet = False
elif alphaPairsChanged == 0:
entireSet = True
print("iteration number: %d" % iter)
return oS
def calcWs(alphas, dataArr, classLabels):
"""
基於alpha計算w值
:param alphas: 拉格朗日乘子
:param dataArr: 數據集
:param classLabels:
:return:
wc 迴歸係數
"""
X = np.mat(dataArr)
labelMat = np.mat(classLabels).transpose()
m, n = np.shape(X)
w = np.zeros((n, 1))
for i in range(m):
w += np.multiply(alphas[i] * labelMat[i], X[i, :].T)
return w
def plotfig_SVM(xArr, yArr, ws, b, alphas):
"""
參考地址：
http://blog.csdn.net/maoersong/article/details/24315633
http://www.cnblogs.com/JustForCS/p/5283489.html
http://blog.csdn.net/kkxgx/article/details/6951959
"""
xMat = np.mat(xArr)
yMat = np.mat(yArr)
# b原來是矩陣，先轉爲數組類型後其數組大小爲（1,1），因此後面加[0]，變爲(1,)
b = np.array(b)[0]
fig = plt.figure()
ax = fig.add_subplot(111)
# 注意flatten的用法
ax.scatter(xMat[:, 0].flatten().A[0], xMat[:, 1].flatten().A[0])
# x最大值，最小值根據原數據集dataArr[:, 0]的大小而定
x = np.arange(-1.0, 10.0, 0.1)
# 根據x.w + b = 0 獲得，其式子展開爲w0.x1 + w1.x2 + b = 0, x2就是y值
y = (-b-ws[0, 0]*x)/ws[1, 0]
ax.plot(x, y)
for i in range(np.shape(yMat[0, :])[1]):
if yMat[0, i] > 0:
ax.plot(xMat[i, 0], xMat[i, 1], 'cx')
else:
ax.plot(xMat[i, 0], xMat[i, 1], 'kp')
# 找到支持向量，並在圖中標紅
for i in range(70):
if alphas[i] > 0.0:
ax.plot(xMat[i, 0], xMat[i, 1], 'ro')
plt.show()
def predict(data, oS):
r = oS.b
for i in range(oS.m):
r += oS.alphas[i] * oS.labelMat[i] * data * oS.X[i, :].T
return 1 if r > 0 else -1
def score(X_test, y_test, oS):
right_count = 0
for i in range(len(X_test)):
result = predict(X_test[i], oS)
if result == y_test[i]:
right_count += 1
return right_count / len(X_test)
if __name__ == "__main__":
# 獲取特徵和目標變量
dataArr, labelArr = loadDataSet('testSetRBF2.txt')
X_train, X_test, y_train, y_test = train_test_split(dataArr, labelArr, test_size=0.3, random_state=4)
# print labelArr
# b是常量值， alphas是拉格朗日乘子
# 0.6 0.001 40 0.5666666666666667
oS = smoP(X_train, y_train, 0.6, 0.001, 200, kTup=('rbf', 10))
b = oS.b
alphas = oS.alphas
svInd = np.nonzero(alphas.A > 0)[0]
sVs = np.mat(X_train)[svInd]
labelSV = np.mat(y_train).transpose()[svInd]
print('/n/n/n')
print('b=', b)
print('alphas[alphas>0]=', alphas[alphas > 0])
print('shape(alphas[alphas > 0])=', np.shape(alphas[alphas > 0]))
for i in range(70):
if alphas[i] > 0:
print(dataArr[i], labelArr[i])
# 畫圖
ws = calcWs(alphas, X_train, y_train)
plotfig_SVM(X_train, y_train, ws, b, alphas)
# print(score(X_test, y_test, oS))
datMat = np.mat(X_test)
lableMat = np.mat(y_test).transpose()
m, n = np.shape(datMat)
right_count = 0
for i in range(m):
kernelEval = kernelTrans(sVs, datMat[i, :], ('rbf', 10))
pre = kernelEval.T * np.multiply(labelSV, alphas[svInd]) + b
sign = 1 if pre > 0 else -1
if sign == labelArr[i]:
right_count += 1
print(right_count / len(X_test))
emmmmmm。。。可能看起來仍是一臉懵逼？？？不要緊，再看看這一份，兩份一塊兒看能夠更好的理解哦
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import random
from sklearn.model_selection import train_test_split
from sklearn.svm import SVC
"""
Author:
ITryagain
Modify:
2019-04-09
"""
def loadDataSet(fileName):
"""
讀取數據
Parameters:
fileName - 文件名
Returns:
dataMat - 數據矩陣
labelMat - 數據標籤
"""
dataMat = []
labelMat = []
fr = open(fileName)
for line in fr.readlines(): # 逐行讀取，濾除空格等
lineArr = line.strip().split('\t')
dataMat.append([float(lineArr[0]), float(lineArr[1])]) # 添加數據
labelMat.append(float(lineArr[2])) # 添加標籤
return np.array(dataMat), np.array(labelMat)
class SVM:
def __init__(self, max_iter=100, kernel='linear'):
self.max_iter = max_iter
self._kernel = kernel
def init_args(self, features, labels):
self.m, self.n = features.shape
self.X = features
self.Y = labels
self.b = 0.0
# 將Ei保存在一個列表裏
self.alpha = np.ones(self.m)
self.E = [self._E(i) for i in range(self.m)]
# 鬆弛變量
self.C = 1.0
def _KKT(self, i):
y_g = self._g(i) * self.Y[i]
if self.alpha[i] == 0:
return y_g >= 1
elif 0 < self.alpha[i] < self.C:
return y_g == 1
else:
return y_g <= 1
# g(x)預測值，輸入xi（X[i]）
def _g(self, i):
r = self.b
for j in range(self.m):
r += self.alpha[j] * self.Y[j] * self.kernel(self.X[i], self.X[j])
return r
# 核函數
def kernel(self, x1, x2):
if self._kernel == 'linear':
return sum([x1[k] * x2[k] for k in range(self.n)])
elif self._kernel == 'poly':
return (sum([x1[k] * x2[k] for k in range(self.n)]) + 1) ** 2
return 0
# E（x）爲g(x)對輸入x的預測值和y的差
def _E(self, i):
return self._g(i) - self.Y[i]
def _init_alpha(self):
# 外層循環首先遍歷全部知足0<a<C的樣本點，檢驗是否知足KKT
index_list = [i for i in range(self.m) if 0 < self.alpha[i] < self.C]
# 不然遍歷整個訓練集
non_satisfy_list = [i for i in range(self.m) if i not in index_list]
index_list.extend(non_satisfy_list)
for i in index_list:
if self._KKT(i):
continue
E1 = self.E[i]
# 若是E2是+，選擇最小的；若是E2是負的，選擇最大的
if E1 >= 0:
j = min(range(self.m), key=lambda x: self.E[x])
else:
j = max(range(self.m), key=lambda x: self.E[x])
return i, j
def _compare(self, _alpha, L, H):
if _alpha > H:
return H
elif _alpha < L:
return L
else:
return _alpha
def fit(self, features, labels):
self.init_args(features, labels)
for t in range(self.max_iter):
# train
i1, i2 = self._init_alpha()
# 邊界
if self.Y[i1] == self.Y[i2]:
L = max(0, self.alpha[i1] + self.alpha[i2] - self.C)
H = min(self.C, self.alpha[i1] + self.alpha[i2])
else:
L = max(0, self.alpha[i2] - self.alpha[i1])
H = min(self.C, self.C + self.alpha[i2] - self.alpha[i1])
E1 = self.E[i1]
E2 = self.E[i2]
# eta=K11+K22-2K12
eta = self.kernel(self.X[i1], self.X[i1]) + self.kernel(self.X[i2], self.X[i2]) - 2 * self.kernel(
self.X[i1], self.X[i2])
if eta <= 0:
# print('eta <= 0')
continue
alpha2_new_unc = self.alpha[i2] + self.Y[i2] * (E1 - E2) / eta # 此處有修改，根據書上應該是E1 - E2，書上130-131頁
alpha2_new = self._compare(alpha2_new_unc, L, H)
alpha1_new = self.alpha[i1] + self.Y[i1] * self.Y[i2] * (self.alpha[i2] - alpha2_new)
b1_new = -E1 - self.Y[i1] * self.kernel(self.X[i1], self.X[i1]) * (alpha1_new - self.alpha[i1]) - self.Y[
i2] * self.kernel(self.X[i2], self.X[i1]) * (alpha2_new - self.alpha[i2]) + self.b
b2_new = -E2 - self.Y[i1] * self.kernel(self.X[i1], self.X[i2]) * (alpha1_new - self.alpha[i1]) - self.Y[
i2] * self.kernel(self.X[i2], self.X[i2]) * (alpha2_new - self.alpha[i2]) + self.b
if 0 < alpha1_new < self.C:
b_new = b1_new
elif 0 < alpha2_new < self.C:
b_new = b2_new
else:
# 選擇中點
b_new = (b1_new + b2_new) / 2
# 更新參數
self.alpha[i1] = alpha1_new
self.alpha[i2] = alpha2_new
self.b = b_new
self.E[i1] = self._E(i1)
self.E[i2] = self._E(i2)
return 'train done!'
def predict(self, data):
r = self.b
for i in range(self.m):
r += self.alpha[i] * self.Y[i] * self.kernel(data, self.X[i])
return 1 if r > 0 else -1
def score(self, X_test, y_test):
right_count = 0
for i in range(len(X_test)):
result = self.predict(X_test[i])
if result == y_test[i]:
right_count += 1
return right_count / len(X_test)
def _weight(self):
# linear model
yx = self.Y.reshape(-1, 1) * self.X
self.w = np.dot(yx.T, self.alpha)
return self.w
if __name__ == '__main__':
X, y = loadDataSet('testSetRBF2.txt')
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=4)
svm = SVM(max_iter=200, kernel='poly')
svm.fit(X, y)
print('the score = {}'.format(svm.score(X_test, y_test)))
clf = SVC(gamma='auto')
clf.fit(X_train, y_train)
print(clf.score(X_test, y_test))
# x_ponits = np.arange(4, 8)
# y_ = -(lr_clf.weights[1] * x_ponits + lr_clf.weights[0]) / lr_clf.weights[2]
# plt.plot(x_ponits, y_)
#
# # lr_clf.show_graph()
data_plus = []
data_minus = []
for i in range(len(X)):
if y[i] > 0:
data_plus.append(X[i])
else:
data_minus.append(X[i])
data_plus_np = np.array(data_plus) # 轉換爲numpy矩陣
data_minus_np = np.array(data_minus) # 轉換爲numpy矩陣
plt.scatter(np.transpose(data_plus_np)[0], np.transpose(data_plus_np)[1]) # 正樣本散點圖
plt.scatter(np.transpose(data_minus_np)[0], np.transpose(data_minus_np)[1]) # 負樣本散點圖
plt.show()
輸出
the score = 0.7 0.8333333333333334
小節
推了三四天的SVM。。。感受本身快死了，若是不是手頭有個比賽須要修改SVM，打死我都不會去推這公式，下次仍是從代碼入手好了，數學很差的痛苦orz。。。。
- 1. 機器學習 支持向量機（SVM） 從理論到放棄，從代碼到理解
- 2. 支持向量機（SVM）從入門到放棄再到掌握
- 3. 機器學習：支持向量機(svm)
- 4. 機器學習：支持向量機(SVM)
- 5. 機器學習：支持向量機SVM
- 6. 機器學習——支持向量機（SVM）
- 7. 支持向量機SVM-----機器學習
- 8. 【機器學習】支持向量機(SVM)
- 9. 機器學習----SVM支持向量機
- 10. 機器學習-支持向量機SVM
- 更多相關文章...
- • Git 分支管理 - Git 教程
- • CAP理論是什麼？ - NoSQL教程
- • 漫談MySQL的鎖機制
- • Docker 清理命令
-
每一个你不满意的现在，都有一个你没有努力的曾经。