ML07-SVM分类

一、SVM模型

1. SVM功能体验

首先通过一个例子来了解SVM的作用；不用关注该例子的代码，直接观察图示效果即可。

数据集显示

% matplotlib inline
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
import sklearn as sk

data = pd.read_csv('train_data.csv')
# feature1  feature2  label
# print(data)
print(data.shape)

# 可视化数据集
figure=plt.figure(figsize=(6,4))
ax=figure.add_axes([0.1,0.1,0.8,0.8],label='数据集图示')
ax.scatter(data['feature1'],data['feature2'],marker='.',color=(1,0,0,1))
plt.show()

(100, 3)

没有分类显示的数据集

数据集分类显示

% matplotlib inline
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
import sklearn as sk

data = pd.read_csv('train_data.csv')
# feature1  feature2  label
# print(data)

#df.sort_values(by="sales" , ascending=False)  #按照某列排序
# 按照标签分类
# print(data['label'].isin([1.0]))
label_1=data[data['label'].isin([1.0])]
label_2=data[data['label'].isin([-1.0])]

# 可视化数据集
figure=plt.figure(figsize=(8,4))
ax=figure.add_axes([0.1,0.1,0.8,0.8],label='数据集图示')
ax.scatter(label_1['feature1'],label_1['feature2'],marker='.',color=(1,0,0,1), label='A类')
ax.scatter(label_2['feature1'],label_2['feature2'],marker='x',color=(0,0,1,1), label='B类')

plt.legend()
plt.show()

数据分类显示结果

分类过程

% matplotlib inline
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
import sklearn as sk
from sklearn import svm

#加载数据集
data = pd.read_csv('train_data.csv')
# feature1  feature2  label
X=data[['feature1','feature2']]
y=data['label']
# 核：线性，
clf = svm.SVC(kernel='rbf', C=0.8, gamma=20)
# 构造gram的计算方式
#gram = np.dot(X, X.T)
# 训练
clf.fit(X, y)

# 预测
predict=clf.predict(X)
print(predict)

print('-------------')
print('SVM Score:',clf.score(X, y))
print('-------------')
# 计算预测值与真实值的相同的个数
count=(predict==y).sum()
print('识别正确数：',count,'总的训练样本数:',y.shape[0])
print('识别正确率：','%5.2f%%'%(100.0*count/y.shape[0]))

[-1.  1. -1.  1. -1. -1. -1.  1. -1.  1.  1. -1.  1.  1.  1. -1.  1. -1.
  1.  1. -1. -1.  1.  1. -1. -1.  1. -1.  1.  1. -1. -1. -1.  1.  1.  1.
 -1. -1. -1.  1.  1. -1. -1. -1.  1.  1. -1.  1. -1.  1.  1.  1.  1. -1.
  1.  1.  1. -1. -1. -1.  1.  1.  1. -1.  1. -1. -1. -1.  1. -1. -1.  1.
 -1. -1.  1.  1. -1.  1. -1. -1.  1. -1. -1.  1. -1. -1.  1. -1. -1. -1.
 -1. -1. -1. -1.  1. -1. -1. -1. -1. -1.]
-------------
SVM Score: 1.0
-------------
识别正确数： 100 总的训练样本数: 100
识别正确率： 100.00%

可视化分类效果

% matplotlib inline
import matplotlib as mpl
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
import sklearn as sk
from sklearn import svm


#加载数据集
data = pd.read_csv('train_data.csv')
# feature1  feature2  label
X=data[['feature1','feature2']]
x1=X['feature1']
x2=X['feature2']
y=data['label']
# 生成网格
x1_min, x1_max = x1.min(), x1.max()  # 第0列的范围
x2_min, x2_max = x2.min(), x2.max()  # 第1列的范围

x1_min=x1_min-0.25
x1_max=x1_max+0.25

x2_min=x2_min-0.25
x2_max=x2_max+0.25

x1, x2 = np.mgrid[x1_min:x1_max:200j, x2_min:x2_max:200j]  # 生成网格采样点
grid_plane = np.stack((x1.flat, x2.flat), axis=1) 
# ----------svm训练开始-----------
# 核：线性，
clf = svm.SVC(kernel='rbf', C=0.8, gamma=20)
# 构造gram的计算方式
#gram = np.dot(X, X.T)
# 训练
clf.fit(X, y)

# 预测
predict=clf.predict(X)

# ----------svm训练结束-----------
z=clf.predict(grid_plane)

# 控制区域颜色
colors = mpl.colors.ListedColormap([(0.2,0.8,0.8,1), (0.8,0.6,0.5,1)])
# 样本数据分开
label_1=data[data['label'].isin([1.0])]
label_2=data[data['label'].isin([-1.0])]

figure=plt.figure(figsize=(8,4))
ax=figure.add_axes([0.1,0.1,0.8,0.8],label='数据集图示')
ax.set_xbound(x1_min,x1_max)
ax.set_ybound(x2_min,x2_max)

ax.pcolormesh(x1, x2,z.reshape(x1.shape),cmap=colors)

ax.scatter(label_1['feature1'],label_1['feature2'],marker='.',color=(1,0,0,1), label='A类')
ax.scatter(label_2['feature1'],label_2['feature2'],marker='x',color=(0,0,1,1), label='B类')

plt.legend()
plt.show()

SVM分类可视化效果

从上面图示看得出来，SVM分类效果非常好，实际上，到目前为止在学界一直公认SVM是最好的分类器。
SVM是一个二元分类算法，线性分类和非线性分类都支持。经过演进，现在也可以支持多元分类，同时经过扩展，也能应用于回归问题。

2. SVM模型

按照机器学习的一般规律，一个智能算法的模型我们分成如下几个部分考虑：
|-决策模型
|-损失函数模型
|-损失函数最小值求解算法模型

一般损失函数的最小值求解算法会有多种，我们再掌握的时候，可以根据个人的工作需要学习。

决策模型

线性分类平面： $xW+b=0$
决策函数： $y=sign(xW+b)$
|- -1
|- 1

损失函数

|-找到一个 $W$ ，使得 $\dfrac{1}{||W||_2}$ 最大，同时使得： $y_i(\vec{x_i} W +b) \ge 1$
|- $L(W,b,\alpha) = \frac{1}{2}|| W || _ 2 ^ 2 - \sum\limits_{i=1}^{m}\alpha_i[y_i(x_i W + b) - 1] \; 满足\alpha_i \geq 0$
其中 $\alpha_i$ 就是拉格朗日乘数。

或者

|- $L(W,b,\xi,\alpha,\mu) = \frac{1}{2}||W||_2^2 +C\sum\limits_{i=1}^{m}\xi_i - \sum\limits_{i=1}^{m}\alpha_i[y_i(x_i W + b) - 1 + \xi_i] - \sum\limits_{i=1}^{m}\mu_i\xi_i$
|-其中 $\alpha_i \geq 0, \mu_i \geq 0$ 是拉格朗日乘数系数。

求解算法

SMO算法

二、SVM的数学基础

1. 从逻辑回归分类开始谈起

下面我们从鸢尾花分类代码开始谈起。

% matplotlib inline
from sklearn.datasets import load_iris
from sklearn.linear_model import LogisticRegression
import matplotlib as  mpl 
import matplotlib.pyplot as plt
import numpy as np

# 鸢尾花的逻辑回归（为了可视化，我们只取两个特征）
X, y = load_iris(return_X_y=True)
X=X[0:100,0:2]
y=y[50:150]

clf = LogisticRegression(random_state=0, solver='lbfgs').fit(X, y)
#clf = LogisticRegression().fit(X, y)
# 预测结果
r=clf.predict(X)

# 权重系数
coef=clf.coef_ 
intercept=clf.intercept_
print('权重系数：',coef)
print('权重截距：',intercept)
score=clf.score(X, y)
print(score)

# 可视化
x1=X[:,0]
x2=X[:,1]

# 生成网格
x1_min, x1_max = x1.min(), x1.max()  # 第0列的范围
x2_min, x2_max = x2.min(), x2.max()  # 第1列的范围

x1_min=x1_min-0.25
x1_max=x1_max+0.25

x2_min=x2_min-0.25
x2_max=x2_max+0.25

x1, x2 = np.mgrid[x1_min:x1_max:200j, x2_min:x2_max:200j]  # 生成网格采样点
grid_plane = np.stack((x1.flat, x2.flat), axis=1) 

z=clf.predict(grid_plane)
# 控制区域颜色
colors = mpl.colors.ListedColormap([(0.2,0.8,0.8,1), (0.8,0.6,0.5,1)])
# 样本数据分开
label_1=y[0:50]
label_2=y[50:100]

figure=plt.figure(figsize=(8,4))
ax=figure.add_axes([0.1,0.1,0.8,0.8],label='数据集图示')
ax.set_xbound(x1_min,x1_max)
ax.set_ybound(x2_min,x2_max)

ax.pcolormesh(x1, x2,z.reshape(x1.shape),cmap=colors)

ax.scatter(X[0:50:,0],X[0:50,1],marker='.',color=(1,0,0,1), label='A类')
ax.scatter(X[50:100:,0],X[50:100,1],marker='x',color=(0,0,1,1), label='B类')

plt.legend()
plt.show()

权重系数： [[ 3.09621805 -3.01680617]]
权重截距： [-7.42914365]
1.0

逻辑回归训练的可视化结果

逻辑回归训练出来的参数为:
|- $W=[[ 3.09621805 -3.01680617]]$
|- $b=[-7.42914365]$

写成函数表达式就是：
|- $y = w_1 x_1 + w_2 x_2 + b$

% matplotlib inline
from sklearn.datasets import load_iris
from sklearn.linear_model import LogisticRegression
import matplotlib as  mpl 
import matplotlib.pyplot as plt
import numpy as np
from sklearn.metrics.pairwise import sigmoid_kernel


# 鸢尾花的逻辑回归（为了可视化，我们只取两个特征）
X, y = load_iris(return_X_y=True)
X=X[0:100,0:2]
y=y[50:150]

clf = LogisticRegression(random_state=0, solver='lbfgs').fit(X, y)
#clf = LogisticRegression().fit(X, y)
# 预测结果
r=clf.predict(X)

# 权重系数
coef=clf.coef_ 
intercept=clf.intercept_
print('权重系数：',coef)
print('权重截距：',intercept)
score=clf.score(X, y)
print(score)

# 可视化
x1=X[:,0]
x2=X[:,1]

# 生成网格
x1_min, x1_max = x1.min(), x1.max()  # 第0列的范围
x2_min, x2_max = x2.min(), x2.max()  # 第1列的范围

x1_min=x1_min-0.25
x1_max=x1_max+0.25

x2_min=x2_min-0.8
x2_max=x2_max+0.8

x1, x2 = np.mgrid[x1_min:x1_max:200j, x2_min:x2_max:200j]  # 生成网格采样点
grid_plane = np.stack((x1.flat, x2.flat), axis=1) 

z=clf.predict(grid_plane)
# 控制区域颜色
colors = mpl.colors.ListedColormap([(0.2,0.8,0.8,1), (0.8,0.6,0.5,1)])
# 样本数据分开
label_1=y[0:50]
label_2=y[50:100]

figure=plt.figure(figsize=(8,6))
ax=figure.add_axes([0.1,0.1,0.8,0.8],label='数据集图示')
ax.set_xbound(x1_min,x1_max)
ax.set_ybound(x2_min,x2_max)

ax.pcolormesh(x1, x2,z.reshape(x1.shape),cmap=colors)

ax.scatter(X[0:50:,0],X[0:50,1],marker='.',color=(1,1,0,1), label='A类样本')
ax.scatter(X[50:100:,0],X[50:100,1],marker='x',color=(0,0,1,1), label='B类样本')

# 按照权重系数绘制一条线
ss=lambda x:  (coef[0,0]*x+intercept)/-coef[0,1]
xx=np.linspace(x1_min,x1_max,200)
yy=ss(xx)
ax.plot(xx,yy,color=(1,0,0,1),linewidth=5,label='逻辑回归线')

plt.legend()
plt.show()

权重系数： [[ 3.09621805 -3.01680617]]
权重截距： [-7.42914365]
1.0

逻辑回归分类线可视化

大家知道逻辑回归的决策函数是：
|- $S( x ) = \dfrac{1}{1-exp^{ (- x W) }}$

从上面图示看出来，一旦加权系数确定后，后面的sigmoid的函数仅仅使用用于输出的数据转换，用来做概率决策使用，不影响分类本身的线性性。（这个理解有点抽象）

根据我们代码，分类线（平面或者超平面）可以表示为：
|- $xW + b = 0$

2. 逻辑回归分类的问题

从上图可以看见，分类面可以多个（如下图），但哪一个是最好的呢？

% matplotlib inline
from sklearn.datasets import load_iris
from sklearn.linear_model import LogisticRegression
import matplotlib as  mpl 
import matplotlib.pyplot as plt
import numpy as np
from sklearn.metrics.pairwise import sigmoid_kernel


# 鸢尾花的逻辑回归（为了可视化，我们只取两个特征）
X, y = load_iris(return_X_y=True)
X=X[0:100,[0,2]]
y=y[50:150]

clf = LogisticRegression(random_state=0, solver='lbfgs').fit(X, y)
#clf = LogisticRegression().fit(X, y)
# 预测结果
r=clf.predict(X)

# 权重系数
coef=clf.coef_ 
intercept=clf.intercept_
print('权重系数：',coef)
print('权重截距：',intercept)
score=clf.score(X, y)
print(score)

# 可视化
x1=X[:,0]
x2=X[:,1]

# 生成网格
x1_min, x1_max = x1.min(), x1.max()  # 第0列的范围
x2_min, x2_max = x2.min(), x2.max()  # 第1列的范围

x1_min=x1_min-0.25
x1_max=x1_max+0.25

x2_min=x2_min-0.8
x2_max=x2_max+0.8

x1, x2 = np.mgrid[x1_min:x1_max:200j, x2_min:x2_max:200j]  # 生成网格采样点
grid_plane = np.stack((x1.flat, x2.flat), axis=1) 

z=clf.predict(grid_plane)
# 控制区域颜色
colors = mpl.colors.ListedColormap([(0.2,0.8,0.8,1), (0.8,0.6,0.5,1)])
# 样本数据分开
label_1=y[0:50]
label_2=y[50:100]

figure=plt.figure(figsize=(8,6))
ax=figure.add_axes([0.1,0.1,0.8,0.8],label='数据集图示')
ax.set_xbound(x1_min,x1_max)
ax.set_ybound(x2_min,x2_max)

ax.pcolormesh(x1, x2,z.reshape(x1.shape),cmap=colors)

ax.scatter(X[0:50:,0],X[0:50,1],marker='.',color=(1,1,0,1), label='A类样本')
ax.scatter(X[50:100:,0],X[50:100,1],marker='x',color=(0,0,1,1), label='B类样本')

# 按照权重系数绘制一条线
ss=lambda x:  (coef[0,0]*x+intercept)/-coef[0,1]
xx=np.linspace(x1_min,x1_max,200)
yy=ss(xx)
ax.plot(xx,yy,color=(0,0,1,1),linewidth=2,label='逻辑回归线')

tt=lambda x:  (coef[0,0]*x+intercept+2.35)/-coef[0,1]
ax.plot(xx,tt(xx),color=(1,0,1,1),linewidth=2,label='支撑平面')

uu=lambda x:  (coef[0,0]*x+intercept-0.75)/-coef[0,1]
ax.plot(xx,uu(xx),color=(1,0,1,1),linewidth=2,label='支撑平面')


vv=lambda x:  (coef[0,0]*x+intercept+0.8)/-coef[0,1]
ax.plot(xx,vv(xx),color=(1,0,0,1),linewidth=5,label='中间平面')

plt.legend()
plt.show()

权重系数： [[0.36638806 2.84169262]]
权重截距： [-9.65298039]
1.0

支撑点可视化

直观上，觉得分类平面在中间可能分类效果更加好。
可以扩展成更加科学的描述：
|-正面描述：找出一个分类平面，使得所有正确分类的样本点到分类平面距离最大；
|-反面描述：找出一个分类平面，使得所有错误分类的样本点到分类平面距离最小。

3. 分类效果的评判标准

分类定义
首先做一个简单的规范定义（或者假设），我们假设分类平面表示为 $\vec{x}W+b=0$ ，在平面的上方假设是1，则在平面下方就为-1。
点到平面（或者直线）的垂直距离表示
假设直线 $w_1 x + w_2 y + b =0$ 为 $\vec{x}W + b = 0$ ，点 $\vec{x_0}=(x_0,y_0)$ 到直线 $\vec{x}W + b = 0$ 的垂直距离为：

|- $d = \dfrac{ | w_1 x_0 + w_2 y_0 + b | }{ \sqrt{ { w_1 }^2 + { w_2 }^2 } } = \dfrac{ | w_1 x_0 + w_2 y_0 + b | }{ || W || _ 2}$
其中 $|| W || _ 2$ 表示二阶范数。

参考百度文库一个关于点到直线的距离证明推导：

https://wenku.baidu.com/view/b123b02291c69ec3d5bbfd0a79563c1ec5dad794.html
函数距离与集合距离
我们把点 $\vec{x}=(x_0,y_0)$ 到直线 $\vec{x} W + b = 0$ 的垂直距离称为几何距离。
我们把 $| w_1 x_0 + w_2 y_0 + b |$ 称为函数距离。
支撑向量
为了让所有样本点到分类平面的距离最大，这样样本被分类平面可分，还与分类平面保持一定距离，这样分类平面就是唯一的。我们可以换成另外一种思维：
找到平行分类平面，且与样本集最近的平面，这个平面对两类样本来说就有两个，这两个平行平面称为支撑向量。

由于所有样本点到分类平面的距离可以成比例放大缩小，所以我们可以假设支撑向量到分类平面向量的函数距离是1，则几何距离就是：
|- $\dfrac{1}{||W||_2}$
支撑向量之间的几何距离是： $\dfrac{2}{|| W || _ 2}$
SVM的优化模型
找到支撑向量，使得其他样本点在支撑向量的两边。
描述成数学方式就是：
|-找到 $W,b$ 与至少一个样本x，使得 $\dfrac{ | \vec{x} W + b | }{ ||W||_2 }=\dfrac{ y(\vec{x} W + b) }{ ||W||_2 }$ 最大,同时使得对每个样本点满足： $y_i(\vec{x_i} W +b) \ge y(\vec{x} W + b)$
其中 $y=1$ 或者 $-1$ ，正负号与 $\vec{x} W + b$ 相同。 $i$ 是第 $i$ 个样本。

未来简化问题，可以固定上式的分子，则可以把支撑向量到分类平面的函数距离假设为1，则上述公司可以表达为。
|-找到一个 $W$ ，使得 $\dfrac{1}{||W||_2}$ 最大，同时使得： $y_i(\vec{x_i} W +b) \ge 1$
可以换一个方向表述：
|-找到一个 $W$ ，使得 $||W||_2$ 最小，同时使得： $y_i(\vec{x_i} W +b) \ge 1$ ，或者

|-找到一个 $W$ ，使得 $\dfrac{1}{2}||W||_2 ^2$ 最小（注意添加了一个平方，这是为了防止开方的麻烦），同时使得： $y_i(\vec{x_i} W +b) \ge 1$

三、SVM的求解推导

1. 拉格朗日乘数法

在数学凸优化理论中，拉格朗日乘数法（以数学家约瑟夫·路易斯·拉格朗日命名）是一种求函数极值方法：变量受一个或多个条件所限制的多元函数。
|-这种方法将一个有n个变量与k个约束条件的最优化问题转换为一个有n + k个变量的方程组的极值问题，其变量不受任何约束。
|-这种方法引入了一种新的标量未知数，即拉格朗日乘数：约束方程的梯度（gradient）的线性组合里每个向量的系数。

此方法的证明牵涉到偏微分，全微分或链法，从而找到能让设出的隐函数的微分为零的未知数的值。

2.SVM的拉格朗日乘数最优问题

拉格朗日乘子函数
$L(W,b,\alpha) = \frac{1}{2}|| W || _ 2 ^ 2 - \sum\limits_{i=1}^{m}\alpha_i[y_i(x_i W + b) - 1] \; 满足\alpha_i \geq 0$
其中 $\alpha_i$ 就是拉格朗日乘数。

补充：

如果假设拉格朗日乘数固定取某个值，则上述问题可以使用梯度下降法求解，下面是Tensorflow的代码

% matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import tensorflow as tf
from sklearn import datasets



# 加载数据
data,target= datasets.load_iris(return_X_y=True)
x=data[0:100,[1,3]]     # 取第一个与第三个特征
y=target[0:100]
# 修改y的值为-1 与 1
y[y==0]=-1

y=y.reshape(y.shape[0],1)

#=========================================
# 初始化输入
x_data = tf.placeholder(shape=[None, 2], dtype=tf.float32)
y_target = tf.placeholder(shape=[None, 1], dtype=tf.float32)

# 创建训练变量
W = tf.Variable(tf.random_normal(shape=[2, 1]))
b  = tf.Variable(tf.random_normal(shape=[1, 1]))

# 定义线性模型
o = tf.add(tf.matmul(x_data, W), b)

# 计算拉格朗日乘数函数（损失函数）
# 假设拉格朗日乘数为0.01
alpha = tf.constant([0.008])

# l2 = tf.reduce_sum(tf.square(W))/2.0
l2 = tf.reduce_sum(tf.matmul(W, tf.transpose(W)))/2.0
# 下面取最大值是只计算争取分类的样本
o_sum = tf.reduce_sum(tf.maximum(0.0, tf.subtract(1., tf.multiply(o, y_target))))
loss = tf.add(o_sum, tf.multiply(alpha, l2))

optimizer = tf.train.GradientDescentOptimizer(0.01)
trainer = optimizer.minimize(loss)

#=========================================
# 训练
sess = tf.Session()
init = tf.global_variables_initializer()
sess.run(init)
# 训练对次
for i in range(20000):
    sess.run(trainer, feed_dict={x_data: x, y_target: y})

v_w = sess.run(W)
v_b = sess.run(b)
print(v_w)
print(v_b)
predict=sess.run(o,feed_dict={x_data: x, y_target: y})
#print(predict)
#=========================================
# 可视化
figure=plt.figure(figsize=(8,4))
ax=figure.add_axes([0.1,0.1,0.8,0.8],label='数据集图示')
x_min=0
x_max=7.5
y_min=0.5
y_max=6
ax.set_xbound(x_min,x_max)
ax.set_ybound(y_min,y_max)

ax.scatter(x[0:50,0],x[0:50,1],marker='.',color=(1,0,0,1), label='A类样本')
ax.scatter(x[50:100,0],x[50:100,1],marker='x',color=(0,1,0,1), label='B类样本')

# 绘制分类面（直线）
ss=lambda x:  (v_w[0,0]*x+v_b[0,0])/-v_w[1,0]
xx=np.linspace(x_min,x_max,200)
yy=ss(xx)
ax.plot(xx,yy,color=(0,0,1,1),linewidth=2,label='SVM分类线')

plt.legend()
plt.show()

[[-1.0834631]
 [ 3.5795624]]
[[0.3453444]]

拉格朗日对偶模型的一个建议算法结果

尽管我们还没有完全按照一次性求解最优的 $W,b,\alpha$ 的算法求解最小值，但我们假设一个拉格朗日乘数（这样求解问题又编程求解 $W,b$ 两个参数，同时算是函数式凸的可导的，使用梯度下降是比较好的），来按照梯度下降算法求解最小值，得到的效果也是比较靠谱的，至少比逻辑回归靠谱。
可以调整上面代码中的alpha值观察分类界面的变化（观察支撑向量与支撑样本点），也可以调整鸢尾花数据的样本特征观察，并自仔细理解我们的优化模型的思想。

使用代码理解支撑向量

% matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import tensorflow as tf
from sklearn import datasets



# 加载数据
data,target= datasets.load_iris(return_X_y=True)
x=data[0:100,[1,3]]     # 取第一个与第三个特征
y=target[0:100]
# 修改y的值为-1 与 1
y[y==0]=-1

y=y.reshape(y.shape[0],1)

#=========================================
# 初始化输入
x_data = tf.placeholder(shape=[None, 2], dtype=tf.float32)
y_target = tf.placeholder(shape=[None, 1], dtype=tf.float32)

# 创建训练变量
W = tf.Variable(tf.random_normal(shape=[2, 1]))
b  = tf.Variable(tf.random_normal(shape=[1, 1]))

# 定义线性模型
o = tf.add(tf.matmul(x_data, W), b)

# 计算拉格朗日乘数函数（损失函数）
# 假设拉格朗日乘数为0.01
alpha = tf.constant([0.008])

# l2 = tf.reduce_sum(tf.square(W))/2.0
l2 = tf.reduce_sum(tf.matmul(W, tf.transpose(W)))/2.0
# 下面取最大值是只计算争取分类的样本
o_sum = tf.reduce_sum(tf.maximum(0.0, tf.subtract(1., tf.multiply(o, y_target))))
loss = tf.add(o_sum, tf.multiply(alpha, l2))

optimizer = tf.train.GradientDescentOptimizer(0.01)
trainer = optimizer.minimize(loss)

#=========================================
# 训练
sess = tf.Session()
init = tf.global_variables_initializer()
sess.run(init)
# 训练对次
for i in range(20000):
    sess.run(trainer, feed_dict={x_data: x, y_target: y})

v_w = sess.run(W)
v_b = sess.run(b)
print(v_w)
print(v_b)
predict=sess.run(o,feed_dict={x_data: x, y_target: y})
#print(predict)
#=========================================
# 可视化
figure=plt.figure(figsize=(8,4))
ax=figure.add_axes([0.1,0.1,0.8,0.8],label='数据集图示')
x_min=0
x_max=7.5
y_min=0.5
y_max=6
ax.set_xbound(x_min,x_max)
ax.set_ybound(y_min,y_max)

ax.scatter(x[0:50,0],x[0:50,1],marker='.',color=(1,0,0,1), label='A类样本')
ax.scatter(x[50:100,0],x[50:100,1],marker='x',color=(0,1,0,1), label='B类样本')

# 绘制分类面（直线）
ss=lambda x:  (v_w[0,0]*x+v_b[0,0])/-v_w[1,0]
xx=np.linspace(x_min,x_max,200)
yy=ss(xx)
ax.plot(xx,yy,color=(0,0,1,1),linewidth=2,label='SVM分类线')

# 绘制支撑向量 
tt=lambda x:  (v_w[0,0]*x+v_b[0,0]+1)/-v_w[1,0]
xx=np.linspace(x_min,x_max,200)
yy=tt(xx)
ax.plot(xx,yy,color=(0,0,1,1),linewidth=2,label='支撑向量',linestyle=':')

uu=lambda x:  (v_w[0,0]*x+v_b[0,0]-1)/-v_w[1,0]
xx=np.linspace(x_min,x_max,200)
yy=uu(xx)
ax.plot(xx,yy,color=(0,0,1,1),linewidth=2,label='支撑向量',linestyle=':')


plt.legend()
plt.show()

[[-1.658036]
 [ 3.903697]]
[[1.5776467]]

SVM拉格朗日对偶模型表示的支撑向量

SVM的拉格朗日优化模型：
$\underbrace{min}_{W,b}\; \underbrace{max}_{\alpha_i \geq 0} L(W,b,\alpha)$
找到最大化拉格朗日函数的拉格朗日乘数，然后再找到最小化拉格朗日函数权重系数 $W$ 与截距 $b$ 。
根据凸优问题推导，可以证明上面拉格朗日乘子函数是存在最小值的，这就是著名的KKT条件。拉格朗日乘子函数满足KKT条件，所以存在最小值。（对数学痴迷的可以自行取证明）
因为拉格朗日乘子函数存在最小值，所以可以把最小值问题转换为如下最大值优化问题：
$\underbrace{max}_{\alpha_i \geq 0} \;\underbrace{min}_{W,b}\; L(W,b,\alpha)$

先找到 $W,b$ 满足最小，然后再找到拉格朗日乘数让拉格朗日乘子函数最大。

3. 求解W与b

求解 $W$
$\dfrac{\partial L}{\partial W} = 0 \;\Rightarrow W^T = \sum\limits_{i=1}^{m}\alpha_iy_ix_i$
求解 $b$
$\dfrac{\partial L}{\partial b} = 0 \;\Rightarrow \sum\limits_{i=1}^{m}\alpha_iy_i = 0$
因为 $b$ 的偏导数没有参数 $b$ ， $b$ 的解可以是任意常数。但这个公式在后面会使用；

4. 求解拉格朗日乘数

简化 $\underbrace{min}_{W,b}\; L(W,b,\alpha)$

$\psi(\alpha) = \underbrace{min}_{W,b}\; L(W,b,\alpha)$
$\begin{align} \psi(\alpha) & = \dfrac{1}{2}||W||_2^2 - \sum\limits_{i=1}^{m}\alpha_i[y_i(x_iW + b) - 1] \\ & = \dfrac{1}{2}WW^T-\sum\limits_{i=1}^{m}\alpha_iy_ix_iW - \sum\limits_{i=1}^{m}\alpha_iy_ib + \sum\limits_{i=1}^{m}\alpha_i \\ & = \dfrac{1}{2}W\sum\limits_{i=1}^{m}\alpha_iy_ix_i -\sum\limits_{i=1}^{m}\alpha_iy_ix_iW - \sum\limits_{i=1}^{m}\alpha_iy_ib + \sum\limits_{i=1}^{m}\alpha_i \\ & = \dfrac{1}{2}W\sum\limits_{i=1}^{m}\alpha_iy_ix_i - W\sum\limits_{i=1}^{m}\alpha_iy_ix_i - \sum\limits_{i=1}^{m}\alpha_iy_ib + \sum\limits_{i=1}^{m}\alpha_i \\ & = - \dfrac{1}{2}W\sum\limits_{i=1}^{m}\alpha_iy_ix_i - \sum\limits_{i=1}^{m}\alpha_iy_ib + \sum\limits_{i=1}^{m}\alpha_i \\ & = - \dfrac{1}{2}W\sum\limits_{i=1}^{m}\alpha_iy_ix_i - b\sum\limits_{i=1}^{m}\alpha_iy_i + \sum\limits_{i=1}^{m}\alpha_i \\ & = -\dfrac{1}{2}(\sum\limits_{i=1}^{m}\alpha_iy_ix_i)^T(\sum\limits_{i=1}^{m}\alpha_iy_ix_i) - b\sum\limits_{i=1}^{m}\alpha_iy_i + \sum\limits_{i=1}^{m}\alpha_i \\ & = -\dfrac{1}{2}\sum\limits_{i=1}^{m}\alpha_iy_ix_i^T\sum\limits_{i=1}^{m}\alpha_iy_ix_i - b\sum\limits_{i=1}^{m}\alpha_iy_i + \sum\limits_{i=1}^{m}\alpha_i \\ & = -\dfrac{1}{2}\sum\limits_{i=1}^{m}\alpha_iy_ix_i^T\sum\limits_{i=1}^{m}\alpha_iy_ix_i + \sum\limits_{i=1}^{m}\alpha_i \\ & = -\dfrac{1}{2}\sum\limits_{i=1,j=1}^{m}\alpha_iy_ix_i^T\alpha_jy_jx_j + \sum\limits_{i=1}^{m}\alpha_i \\ & = \sum\limits_{i=1}^{m}\alpha_i - \dfrac{1}{2}\sum\limits_{i=1,j=1}^{m}\alpha_i\alpha_jy_iy_jx_i^Tx_j \end{align}$

求解拉格朗日乘数

对上面简化的函数最大值：
$\psi(\alpha) = \sum\limits_{i=1}^{m}\alpha_i - \dfrac{1}{2}\sum\limits_{i=1,j=1}^{m}\alpha_i\alpha_jy_iy_jx_i^Tx_j$ ，并且同时满足： $\sum\limits_{i=1}^{m}\alpha_iy_i = 0$

转换为最小值问题：
$\psi(\alpha) = \dfrac{1}{2}\sum\limits_{i=1,j=1}^{m}\alpha_i\alpha_jy_iy_jx_i^Tx_j - \sum\limits_{i=1}^{m}\alpha_i$ ，并且同时满足： $\sum\limits_{i=1}^{m}\alpha_iy_i = 0$

求解上述最小值，需要用到SMO（Sequential Minimal Optimization）求解算法。

5. 软距离优化模型

SVM分类中的奇异点问题
下面通过一个例子代码来说明SVM在非完全线性可分下的问题（使用鸢尾花其中第二类与第三类是线性不可分的）。

% matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import tensorflow as tf
from sklearn import datasets



# 加载数据
data,target= datasets.load_iris(return_X_y=True)
x=data[50:150,[1,3]]     # 取第一个与第三个特征
y=target[50:150]
# 修改y的值为-1 与 1
y-=1
y[y==0]=-1

y=y.reshape(y.shape[0],1)

#=========================================
# 初始化输入
x_data = tf.placeholder(shape=[None, 2], dtype=tf.float32)
y_target = tf.placeholder(shape=[None, 1], dtype=tf.float32)

# 创建训练变量
W = tf.Variable(tf.random_normal(shape=[2, 1]))
b  = tf.Variable(tf.random_normal(shape=[1, 1]))

# 定义线性模型
o = tf.add(tf.matmul(x_data, W), b)

# 计算拉格朗日乘数函数（损失函数）
# 假设拉格朗日乘数为0.01
alpha = tf.constant([0.008])

# l2 = tf.reduce_sum(tf.square(W))/2.0
l2 = tf.reduce_sum(tf.matmul(W, tf.transpose(W)))/2.0
# 下面取最大值是只计算争取分类的样本
o_sum = tf.reduce_sum(tf.maximum(0.0, tf.subtract(1., tf.multiply(o, y_target))))
loss = tf.add(o_sum, tf.multiply(alpha, l2))

optimizer = tf.train.GradientDescentOptimizer(0.01)
trainer = optimizer.minimize(loss)

#=========================================
# 训练
sess = tf.Session()
init = tf.global_variables_initializer()
sess.run(init)
# 训练对次
for i in range(20000):
    sess.run(trainer, feed_dict={x_data: x, y_target: y})

v_w = sess.run(W)
v_b = sess.run(b)
print(v_w)
print(v_b)
predict=sess.run(o,feed_dict={x_data: x, y_target: y})
#=========================================
# 可视化
figure=plt.figure(figsize=(8,4))
ax=figure.add_axes([0.1,0.1,0.8,0.8],label='数据集图示')
x_min=0
x_max=7.5
y_min=0.5
y_max=6
ax.set_xbound(x_min,x_max)
ax.set_ybound(y_min,y_max)

ax.scatter(x[0:50,0],x[0:50,1],marker='.',color=(1,0,0,1), label='A类样本')
ax.scatter(x[50:100,0],x[50:100,1],marker='x',color=(0,1,0,1), label='B类样本')

# 绘制分类面（直线）
ss=lambda x:  (v_w[0,0]*x+v_b[0,0])/-v_w[1,0]
xx=np.linspace(x_min,x_max,200)
yy=ss(xx)
ax.plot(xx,yy,color=(0,0,1,1),linewidth=2,label='SVM分类线')

plt.legend()
plt.show()

[[-1.8957276]
 [ 8.625281 ]]
[[-9.2072735]]

异常点对分类的影响的可视化说明

注意上面例子中，某些点在分类中属于异常点，这些异常点可能干扰训练过程，影响识别结果。所以为了避免这种情况，在SVM的优化条件中降低要求，提出软距离的概念来降低干扰。原来提出的支撑向量到分类平面的距离为硬距离，硬距离上面假设是1，所谓软距离就是引入一个修正量，使得支撑向量到分类平面的距离不是绝对大于1，而是引入更加宽松的约束。

软距离优化模型
引入软距离来修正SVM优化问题。
|-找到分类平面，使得 $\frac{1}{2}||W||_2^2 +C\sum\limits_{i=1}^{m}\xi_i$ 最小
|-同时满足约束条件：
|-条件一： $y_i(x_i W + b) \geq 1 - \xi_i \;\;(i =1,2,...m)$
|-条件二： $\xi_i \geq 0 \;\;(i =1,2,...m)$

6. 软距离优化模型的优化目标：拉格朗日乘数函数

拉格朗日乘数函数

$L(W,b,\xi,\alpha,\mu) = \frac{1}{2}||W||_2^2 +C\sum\limits_{i=1}^{m}\xi_i - \sum\limits_{i=1}^{m}\alpha_i[y_i(x_i W + b) - 1 + \xi_i] - \sum\limits_{i=1}^{m}\mu_i\xi_i$
其中 $\alpha_i \geq 0, \mu_i \geq 0$ 是拉格朗日乘数系数。

拉格朗日优化目标

SVM优化目标就是求解拉格朗日函数的最小值，可以表述为如下数学描述：
|- $\underbrace{min}_{W,b,\xi}\; \underbrace{max}_{\alpha_i \geq 0, \mu_i \geq 0,} L(W,b,\alpha, \xi,\mu)$

拉格朗日优化目标转换

数学上已经证明这个函数满足凸与可导性质，而且满足KKT条件：
$\begin{align} \nabla_x L(W,b,\xi,\alpha,\mu) &= 0 \\ C\sum\limits_{i=1}^{m}\xi_i &= 0 \\ y_i(x_i W + b) &\geq 1 - \xi_i \;\;(i =1,2,...m) \\ \xi_i &\geq 0 \;\;(i =1,2,...m) \\ \alpha _i, \mu_i &\ge 0 , \ i=1,2,...,n \\ \end{align}$

其中第一项是求极小值，第二项是松弛条件，第三、四项是约束条件，第五项是拉格朗日的乘数条件

所以，可以如下描述拉格朗日优化目标：
$\underbrace{max}_{\alpha_i \geq 0, \mu_i \geq 0,} \; \underbrace{min}_{W,b,\xi}\; L(W,b,\alpha, \xi,\mu)$

7. 软距离优化模型的推导

首先对 $W,b,\xi$ 求偏导：

$W$ 的偏导：
$\dfrac{\partial L}{\partial W} = 0 \;\Rightarrow W^T = \sum\limits_{i=1}^{m}\alpha_iy_ix_i$
$b$ 的偏导：
$\dfrac{\partial L}{\partial b} = 0 \;\Rightarrow \sum\limits_{i=1}^{m}\alpha_iy_i = 0$
$\xi$ 的偏导：
$\dfrac{\partial L}{\partial \xi} = 0 \;\Rightarrow C- \alpha_i - \mu_i = 0$
消去 $W,b,\xi$ 的拉格朗日乘数函数
$\begin{align} L(W,b,\xi,\alpha,\mu) & = \frac{1}{2}||W||_2^2 +C\sum\limits_{i=1}^{m}\xi_i - \sum\limits_{i=1}^{m}\alpha_i[y_i(x_iW + b) - 1 + \xi_i] - \sum\limits_{i=1}^{m}\mu_i\xi_i 　\\ & = \frac{1}{2}||W||_2^2 +( \alpha_i + \mu_i)\sum\limits_{i=1}^{m}\xi_i - \sum\limits_{i=1}^{m}\alpha_i[y_i(x_iW + b) - 1 + \xi_i] - \sum\limits_{i=1}^{m}\mu_i\xi_i 　\\ & = \frac{1}{2}||W||_2^2 + \alpha_i \sum\limits_{i=1}^{m}\xi_i + \mu_i \sum\limits_{i=1}^{m}\xi_i - \sum\limits_{i=1}^{m}\alpha_i[y_i(x_iW + b) - 1 + \xi_i] - \sum\limits_{i=1}^{m}\mu_i\xi_i 　\\ &= \frac{1}{2}||W||_2^2 - \sum\limits_{i=1}^{m}\alpha_i[y_i(x_iW + b) - 1 + \xi_i] + \sum\limits_{i=1}^{m}\alpha_i\xi_i \\ & = \frac{1}{2}||W||_2^2 - \sum\limits_{i=1}^{m}\alpha_i[y_i(x_iW + b) - 1] \\ & = \frac{1}{2}WW^T-\sum\limits_{i=1}^{m}\alpha_iy_ix_iW - \sum\limits_{i=1}^{m}\alpha_iy_ib + \sum\limits_{i=1}^{m}\alpha_i \\ & = \frac{1}{2}W\sum\limits_{i=1}^{m}\alpha_iy_ix_i -\sum\limits_{i=1}^{m}\alpha_iy_ix_iW - \sum\limits_{i=1}^{m}\alpha_iy_ib + \sum\limits_{i=1}^{m}\alpha_i \\ & = \frac{1}{2}W\sum\limits_{i=1}^{m}\alpha_iy_ix_i - W\sum\limits_{i=1}^{m}\alpha_iy_ix_i - \sum\limits_{i=1}^{m}\alpha_iy_ib + \sum\limits_{i=1}^{m}\alpha_i \\ & = - \frac{1}{2}W\sum\limits_{i=1}^{m}\alpha_iy_ix_i - \sum\limits_{i=1}^{m}\alpha_iy_ib + \sum\limits_{i=1}^{m}\alpha_i \\ & = - \frac{1}{2}W\sum\limits_{i=1}^{m}\alpha_iy_ix_i - b\sum\limits_{i=1}^{m}\alpha_iy_i + \sum\limits_{i=1}^{m}\alpha_i \\ & = -\frac{1}{2}(\sum\limits_{i=1}^{m}\alpha_iy_ix_i)^T(\sum\limits_{i=1}^{m}\alpha_iy_ix_i) - b\sum\limits_{i=1}^{m}\alpha_iy_i + \sum\limits_{i=1}^{m}\alpha_i \\ & = -\frac{1}{2}\sum\limits_{i=1}^{m}\alpha_iy_ix_i^T\sum\limits_{i=1}^{m}\alpha_iy_ix_i - b\sum\limits_{i=1}^{m}\alpha_iy_i + \sum\limits_{i=1}^{m}\alpha_i \\ & = -\frac{1}{2}\sum\limits_{i=1}^{m}\alpha_iy_ix_i^T\sum\limits_{i=1}^{m}\alpha_iy_ix_i + \sum\limits_{i=1}^{m}\alpha_i \\ & = -\frac{1}{2}\sum\limits_{i=1,j=1}^{m}\alpha_iy_ix_i^T\alpha_jy_jx_j + \sum\limits_{i=1}^{m}\alpha_i \\ & = \sum\limits_{i=1}^{m}\alpha_i - \frac{1}{2}\sum\limits_{i=1,j=1}^{m}\alpha_i\alpha_jy_iy_jx_i^Tx_j \end{align}$

8. SVM优化目标简化

优化目标
通过上述的简化，SVM优化目标基本上可以描述为：
|- 找到一个 $\alpha$ 使得： $\sum\limits_{i=1}^{m}\alpha_i - \dfrac{1}{2}\sum\limits_{i=1,j=1}^{m}\alpha_i\alpha_jy_iy_jx_i^Tx_j$ 最大
或者
|- $\underbrace{ max }_{\alpha} \sum\limits_{i=1}^{m}\alpha_i - \dfrac{1}{2}\sum\limits_{i=1,j=1}^{m}\alpha_i\alpha_jy_iy_jx_i^Tx_j$
同时满足几个约束：
|- $\sum\limits_{i=1}^{m}\alpha_iy_i = 0$
|- $C- \alpha_i - \mu_i = 0$
|- $\alpha_i \geq 0 \;(i =1,2,...,m)$
|- $\mu_i \geq 0 \;(i =1,2,...,m)$
上面几个约束，实际上与最开始的约束不一样，是可以带入消除的。
消除相关约束的优化目标：
同时把问题转化为求最小值。
$\underbrace{ min }_{\alpha} \frac{1}{2}\sum\limits_{i=1,j=1}^{m}\alpha_i\alpha_jy_iy_jx_i^Tx_j - \sum\limits_{i=1}^{m}\alpha_i$
|-条件一： $\sum\limits_{i=1}^{m}\alpha_iy_i = 0$
|-条件二： $0 \leq \alpha_i \leq C$
求解方法
通过SMO算法来求上式极小化时对应的 $\alpha$ 就可以求出 $W$ 和 $b$ 了。
大家可以注意，没有松弛约束条件前， $\alpha_i$ 只需要大于0，现在需要小于 $C$ 。 $C$ 这个参数会产生什么影响，大家可以从最早提出C的背景去理解。

9. SVM核

线性回归与非线性回归基本思想

把低维非线性特征转换成高位线性回归特征。下面通过一个例子可以说明这种转换思想。

比如二维的非线性曲线函数表达： $f_a(x_1, x_2) = a_0 + a_{1}x_1 + a_{2}x_{2} +a_{3}x_1^{2} + a_{4}x_2^{2} + a_{5}x_{1}x_2$

假设令： $y_1 = x_1, y_2 = x_2, y_3 =x_1^{2}, y_4 = x_2^{2}, y_5 = x_{1}x_2$

则得到5维线性直线函数表达： $f_a(y_1, y_2, y_3, y_4, y_5) = a_0 + a_{1}y_1 + a_{2}y_{2} + a_{3}y_3 + a_{4}y_4 + a_{5}y_5$

这样非线性可分变成线性可分。我们可以记 $(x_1, x_2)$ 到 $(y_1, y_2, y_3, y_4, y_5)$ 的映射为函数:
|- $y = \phi ( x )$ 其中 $x=(x_1, x_2),\; \; y=(y_1, y_2, y_3, y_4, y_5)$

低维到高位的问题

如果在SVM中，引入低维到高位的思想，也可以使得SVM具有非线性可分能力，则我们上述的SVM软距离优化模式问题可以表述为：
|- $\underbrace{ min }_{\alpha} \dfrac{1}{2}\sum\limits_{i=1,j=1}^{m}\alpha_i\alpha_jy_iy_j\phi(x_i) \phi(x_j) - \sum\limits_{i=1}^{m}\alpha_i$
|-且 $\sum\limits_{i=1}^{m}\alpha_iy_i = 0$
|-且 $0 \leq \alpha_i \leq C$

但是，如果低维维数很多，则高维的计算就会变得非常庞大，如果是对个样本进行计算，这个运算量基本上是无法承受之重。
数学总是解决问题的，有数学家提出了核函数的概念；

核函数

假设 $\phi$ 是一个从低维空间 $X$ （欧式空间的子空间）到高维的希尔伯特空间的 $H$ 映射。那么如果存在函数 $K(x,z)$ ，对于任意 $x,z\in X$ 都有：
|- $K(x, z) = \phi(x_i) \phi(x_j)$

我们称 $K(x, z)$ 为核函数。这样就把大量的高维矩阵运算变成低维空间的函数运算（可以当成函数运算的数据预处理）。

核函数是需要满足一定条件的：正定核函数（关于正定的概念，这里不介绍了，大家可以查找相关的数学概念）。找到一个正定核函数是比较难的。

常见核函数

线性核函数： $K(x, z) = x z$

多项式核函数： $K(x, z) = （\gamma x z + r)^d$ 其中 $\gamma, r,d$ 是用户定义参数，需要反复调试。

高斯核函数（Gaussian Kernel）： $K(x, z) = exp(-\gamma||x-z||^2)$ 其中 $\gamma$ 用户自己设置。
在SVM中也称为径向基核函数（Radial Basis Function,RBF），它是非线性分类SVM最主流的核函数：

sigmoid核函数： $K(x, z) = \dfrac{1}{1+exp^{-（\gamma x \bullet z + r）}}=tanh（\gamma x \bullet z + r) = S（\gamma x \bullet z + r)$
也是线性不可分SVM常用的核函数之一

10. SMO算法

SMO算法的基本思想

SMO算法就是每次选择其中2个 $\alpha$ 参数训练，其他当成常量。然后再选择下两个参数训练。假设是第一次训练，则优化可以做如下处理：
|- 记： $K_{ij} = \phi(x_i) \phi(x_j)$
|- $\;\underbrace{ min }_{\alpha_1, \alpha_1} \frac{1}{2}K_{11}\alpha_1^2 + \frac{1}{2}K_{22}\alpha_2^2 +y_1y_2K_{12}\alpha_1 \alpha_2 -(\alpha_1 + \alpha_2) +y_1\alpha_1\sum\limits_{i=3}^{m}y_i\alpha_iK_{i1} + y_2\alpha_2\sum\limits_{i=3}^{m}y_i\alpha_iK_{i2}$
|- 且： $\alpha_1y_1 + \alpha_2y_2 = -\sum\limits_{i=3}^{m}y_i\alpha_i = \varsigma$
|- 且： $0 \leq \alpha_i \leq C \;\; i =1,2$

SMO算法的计算公式

（1）先得到： $\alpha_2^{new,unc} = \alpha_2^{old} + \dfrac{y_2(E_1-E_2)}{K_{11} +K_{22}-2K_{12})}$
（2）然后根据 $\alpha_2$ 与 $\alpha_1$ 的关系 $\alpha_1 = y_1(\varsigma - \alpha_2y_2)$ 得到 $\alpha_1$
（3） $\varsigma=\alpha_1y_1 + \alpha_2y_2 = -\sum\limits_{i=3}^{m}y_i\alpha_i$

11. SVM算法总结

决定优化目标

$\underbrace{min}_{\alpha} \frac{1}{2}\sum\limits_{i=1}^{m}\sum\limits_{j=1}^{m}\alpha_i\alpha_jy_iy_j(x_i \bullet x_j) - \sum\limits_{i=1}^{m} \alpha_i$
且： $\sum\limits_{i=1}^{m}\alpha_iy_i = 0$
且： $\alpha_i \geq 0 \; i=1,2,...m$

求解对偶系数

用SMO算法求出所有样本最小时对应的 $\alpha_i\;\;其中i \in (1,2,\dots,m)$

求解权重系数

根据已知的 $\alpha_i，其中i \in(1,2,\dots,m)$ ，求解 $W$ ： $W = \sum\limits_{i=1}^{m}\alpha_i^{*}y_ix_i$

找到支撑向量

找到支撑向量的办法实际比较简单，只需要 $\alpha_i > 0$ 且满足 $y_i(w^Tx_i + b) - 1 = 0$ 对样的样本点就是支撑向量。

根据 $y_i(x_i W + b) \geq 1 (i =1,2,...m)$ 条件，在拉格朗日对偶互换得到 $\alpha_{i}(y_i(w^Tx_i + b) - 1) \geq 0$

满足 $\alpha_{i}(y_i(w^Tx_i + b) - 1) = 0$ 的样本点就是支撑向量：
|- 情况一： $\alpha_i>0$ ，支撑向量必须满足： $y_i(w^Tx_i + b) - 1) = 0$
|- 情况二： $\alpha_i=0$ ，则对应的样本要么是支撑向量，要么就是被正确分类。

根据支撑向量计算截距

根据支撑向量点，就可以计算出解决：
|-根据支撑向量计算满足的条件： $y_s(\sum\limits_{i=1}^{m}\alpha_iy_ix_i^Tx_s+b) = 1$ ， $s \in$ {支撑向量集}，计算支撑向量对应的截距 $b_s$ ， $b_s = y_s - \sum\limits_{i=1}^{s}\alpha_iy_ix_i^Tx_s$
|-根据所有支撑向量的 $b_s$ 计算出 $b$ ： $b = \dfrac{1}{S}\sum\limits_{i=1}^{S}b_s$ ，其中 $S$ 是支撑向量样本数。

四、SVM的实现与分类应用

重点理解各个参数的来源与作用。

1. 线性可分应用

% matplotlib inline
import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
import sklearn as sk
from sklearn import svm
from sklearn import datasets

#加载数据集
data,target = datasets.load_iris(return_X_y=True)
# feature1  feature2  label
X=data[0:100, [0,2]]
x1=X[:,0]
x2=X[:,1]
y=target[0:100]

# 把标签中为0的替换成-1
y[y==0]=-1
# 生成网格
x_min=4
x_max=7.5
y_min=0.5
y_max=6



g1, g2 = np.mgrid[x_min:x_max:200j, y_min:y_max:200j]  # 生成网格采样点
grid_plane = np.stack((g1.flat, g2.flat), axis=1)     # 生成网格，每个点一个坐标
# ----------svm训练开始-----------
# 核：线性，
clf = svm.SVC(kernel='linear', C=3, gamma=20, probability=True)
# 构造gram的计算方式
#gram = np.dot(X, X.T)
# 训练
clf.fit(X, y)
# ----------svm训练结束-----------


# ----------预测开始-----------
predict=clf.predict(X)
#print(predict)
#print(clf.predict_log_proba(X))
#print(clf.predict_proba(X))
# ----------预测结束-----------

# 几个重要的属性
print('支撑向量索引：\n',clf.support_)
print('支撑向量：\n',clf.support_vectors_)
print('每一类支撑向量个数：\n',clf.n_support_)
print('支撑向量的对偶系数：\n',clf.dual_coef_)
print('权重：\n',clf.coef_)
print('截距：\n',clf.intercept_)
print('0（成功）或者1（失败）：\n',clf.fit_status_)
print('A类概率评估：\n',clf.probA_)
print('B类概率评估：\n',clf.probB_)
print('评估正确率：\n',clf.score(X,y))
# ----------svm训练结束-----------
# 绘制区域网格的
z=clf.predict(grid_plane)
# 控制区域颜色
colors = mpl.colors.ListedColormap([(0.2,0.8,0.8,1), (0.8,0.6,0.5,1)])
# ----------区域网格数据结束-----------

figure=plt.figure(figsize=(8,4))
ax=figure.add_axes([0.1,0.1,0.8,0.8],label='数据集图示')


ax.pcolormesh(g1, g2,z.reshape(g1.shape),cmap=colors)

ax.scatter(x1[0:50],x2[0:50],marker='.',color=(1,0,0,1), label='A类')
ax.scatter(x1[50:100],x2[50:100],marker='x',color=(0,0,1,1), label='B类')

#  绘制支撑向量
ax.scatter(clf.support_vectors_[:,0],clf.support_vectors_[:,1],marker='8',color=(0,1,0,1), label='支撑向量', )

# 支撑向量平面
# 绘制支撑向量 

ax.set_xbound(x_min,x_max)
ax.set_ybound(y_min,y_max)

tt=lambda x:  (clf.coef_[0,0]*x+clf.intercept_+1)/-clf.coef_[0,1]
xx=np.linspace(x_min,x_max,200)
yy=tt(xx)
ax.plot(xx,yy,color=(0,0,1,1),linewidth=2,label='支撑向量',linestyle=':')

uu=lambda x:  (clf.coef_[0,0]*x+clf.intercept_-1)/-clf.coef_[0,1]
xx=np.linspace(x_min,x_max,200)
yy=uu(xx)
ax.plot(xx,yy,color=(0,0,1,1),linewidth=2,label='支撑向量',linestyle=':')

plt.legend()
plt.show()

支撑向量索引：
 [24 44 98]
支撑向量：
 [[4.8 1.9]
 [5.1 1.9]
 [5.1 3. ]]
每一类支撑向量个数：
 [2 1]
支撑向量的对偶系数：
 [[-0.00393623 -1.64866481  1.65260103]]
权重：
 [[1.18086785e-03 1.81786114e+00]]
截距：
 [-4.45972162]
0（成功）或者1（失败）：
 0
A类概率评估：
 [-1.72969438]
B类概率评估：
 [-0.63251932]
评估正确率：
 1.0

线性可分训练结果

其中明显可以从支撑点看见软距离导致的惩罚效果。大家可以修改C观察效果。

其中 $probA,probB$ 可以用于如下分类概率公式 $P(y=1 \| f(x))=\dfrac{1}{1+e^{probAf(x)+probB}}$ 。

dual_coef_是支撑向量机对应的对偶系数。

2. 非线性可分应用

实际上除了线性SVM核，其他都是非线性可分的。

多项式核(poly)
大家可以比较下非线性与线性核的计算速度差异。

% matplotlib inline
import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
import sklearn as sk
from sklearn import svm
from sklearn import datasets

#加载数据集
data,target = datasets.load_iris(return_X_y=True)
# feature1  feature2  label
X=data[50:150, [0,2]]
x1=X[:,0]
x2=X[:,1]
y=target[0:100]

# 把标签中为0的替换成-1
y[y==0]=-1
# 生成网格
x_min=4
x_max=10
y_min=2.5
y_max=7.5



g1, g2 = np.mgrid[x_min:x_max:200j, y_min:y_max:200j]  # 生成网格采样点
grid_plane = np.stack((g1.flat, g2.flat), axis=1)     # 生成网格，每个点一个坐标
# ----------svm训练开始-----------
# 核：线性，
clf = svm.SVC(kernel='poly', C=3, gamma=20, probability=True)
# 构造gram的计算方式
#gram = np.dot(X, X.T)
# 训练
clf.fit(X, y)
# ----------svm训练结束-----------


# ----------预测开始-----------
predict=clf.predict(X)
#print(predict)
#print(clf.predict_log_proba(X))
#print(clf.predict_proba(X))
# ----------预测结束-----------

# 几个重要的属性
print('支撑向量索引：\n',clf.support_)
print('支撑向量：\n',clf.support_vectors_)
print('每一类支撑向量个数：\n',clf.n_support_)
print('支撑向量的对偶系数：\n',clf.dual_coef_)
#  print('权重：\n',clf.coef_)     #对非线性，不存在
print('截距：\n',clf.intercept_)
print('0（成功）或者1（失败）：\n',clf.fit_status_)
print('A类概率评估：\n',clf.probA_)
print('B类概率评估：\n',clf.probB_)
print('评估正确率：\n',clf.score(X,y))
# ----------svm训练结束-----------
# 绘制区域网格的
z=clf.predict(grid_plane)
# 控制区域颜色
colors = mpl.colors.ListedColormap([(0.2,0.8,0.8,1), (0.8,0.6,0.5,1)])
# ----------区域网格数据结束-----------

figure=plt.figure(figsize=(8,8))
ax=figure.add_axes([0.1,0.1,0.8,0.8],label='数据集图示')


ax.pcolormesh(g1, g2,z.reshape(g1.shape),cmap=colors)

ax.scatter(x1[0:50],x2[0:50],marker='.',color=(1,0,0,1), label='A类')
ax.scatter(x1[50:100],x2[50:100],marker='x',color=(0,0,1,1), label='B类')

#  绘制支撑向量
ax.scatter(clf.support_vectors_[:clf.n_support_[0],0],clf.support_vectors_[:clf.n_support_[0],1],marker='8',color=(0,1,0,1), label='A类支撑向量', )
ax.scatter(clf.support_vectors_[clf.n_support_[0]:,0],clf.support_vectors_[clf.n_support_[0]:,1],marker='8',color=(1,0,1,1), label='B类支撑向量', )
# 支撑向量平面
# 绘制支撑向量 

ax.set_xbound(x_min,x_max)
ax.set_ybound(y_min,y_max)

plt.legend()
plt.show()

支撑向量索引：
 [ 2  5 22 33 34 56 76 88 91 96]
支撑向量：
 [[6.9 4.9]
 [5.7 4.5]
 [6.3 4.9]
 [6.  5.1]
 [5.4 4.5]
 [4.9 4.5]
 [6.2 4.8]
 [6.  4.8]
 [6.9 5.1]
 [6.3 5. ]]
每一类支撑向量个数：
 [5 5]
支撑向量的对偶系数：
 [[-2.11130059 -3.         -3.         -3.         -0.38922875  2.01937859
   2.81893845  3.          3.          0.6622123 ]]
截距：
 [-4303.16049072]
0（成功）或者1（失败）：
 0
A类概率评估：
 [-0.00178792]
B类概率评估：
 [-0.3176556]
评估正确率：
 0.95

SVM多项式内核

考虑degree与惩罚系数小于1的情况，同时考虑gamma参数值。

% matplotlib inline
import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
import sklearn as sk
from sklearn import svm
from sklearn import datasets

#加载数据集
data,target = datasets.load_iris(return_X_y=True)
# feature1  feature2  label
X=data[50:150, [0,2]]
x1=X[:,0]
x2=X[:,1]
y=target[0:100]

# 把标签中为0的替换成-1
y[y==0]=-1
# 生成网格
x_min=4
x_max=10
y_min=2.5
y_max=7.5

g1, g2 = np.mgrid[x_min:x_max:200j, y_min:y_max:200j]  # 生成网格采样点
grid_plane = np.stack((g1.flat, g2.flat), axis=1)     # 生成网格，每个点一个坐标
# ----------svm训练开始-----------
# 核：线性，
clf = svm.SVC(kernel='poly', C=100, gamma=20, probability=True,degree=3)
# 构造gram的计算方式
#gram = np.dot(X, X.T)
# 训练
clf.fit(X, y)
# ----------svm训练结束-----------


# ----------预测开始-----------
predict=clf.predict(X)
#print(predict)
#print(clf.predict_log_proba(X))
#print(clf.predict_proba(X))
# ----------预测结束-----------

# 几个重要的属性
print('支撑向量索引：\n',clf.support_)
print('支撑向量：\n',clf.support_vectors_)
print('每一类支撑向量个数：\n',clf.n_support_)
print('支撑向量的对偶系数：\n',clf.dual_coef_)
#  print('权重：\n',clf.coef_)     #对非线性，不存在
print('截距：\n',clf.intercept_)
print('0（成功）或者1（失败）：\n',clf.fit_status_)
print('A类概率评估：\n',clf.probA_)
print('B类概率评估：\n',clf.probB_)
print('评估正确率：\n',clf.score(X,y))
# ----------svm训练结束-----------
# 绘制区域网格的
z=clf.predict(grid_plane)
# 控制区域颜色
colors = mpl.colors.ListedColormap([(0.2,0.8,0.8,1), (0.8,0.6,0.5,1)])
# ----------区域网格数据结束-----------

figure=plt.figure(figsize=(8,8))
ax=figure.add_axes([0.1,0.1,0.8,0.8],label='数据集图示')


ax.pcolormesh(g1, g2,z.reshape(g1.shape),cmap=colors)

ax.scatter(x1[0:50],x2[0:50],marker='.',color=(1,0,0,1), label='A类')
ax.scatter(x1[50:100],x2[50:100],marker='x',color=(0,0,1,1), label='B类')

#  绘制支撑向量
ax.scatter(clf.support_vectors_[:clf.n_support_[0],0],clf.support_vectors_[:clf.n_support_[0],1],marker='8',color=(0,1,0,1), label='A类支撑向量', )
ax.scatter(clf.support_vectors_[clf.n_support_[0]:,0],clf.support_vectors_[clf.n_support_[0]:,1],marker='8',color=(1,0,1,1), label='B类支撑向量', )
# 支撑向量平面
# 绘制支撑向量 

ax.set_xbound(x_min,x_max)
ax.set_ybound(y_min,y_max)

plt.legend()
plt.show()

支撑向量索引：
 [ 2  5 22 33 34 56 76 88 91 96]
支撑向量：
 [[6.9 4.9]
 [5.7 4.5]
 [6.3 4.9]
 [6.  5.1]
 [5.4 4.5]
 [4.9 4.5]
 [6.2 4.8]
 [6.  4.8]
 [6.9 5.1]
 [6.3 5. ]]
每一类支撑向量个数：
 [5 5]
支撑向量的对偶系数：
 [[ -70.26955078 -100.         -100.         -100.          -12.60458796
    67.16945263   93.37821998  100.          100.           22.32646612]]
截距：
 [-142944.10942419]
0（成功）或者1（失败）：
 0
A类概率评估：
 [-4.65437252e-05]
B类概率评估：
 [-0.15516355]
评估正确率：
 0.95

degree对多项式的影响效果

高斯核(rbf)

% matplotlib inline
import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
import sklearn as sk
from sklearn import svm
from sklearn import datasets

#加载数据集
data,target = datasets.load_iris(return_X_y=True)
# feature1  feature2  label
X=data[50:150, [0,2]]
x1=X[:,0]
x2=X[:,1]
y=target[0:100]

# 把标签中为0的替换成-1
y[y==0]=-1
# 生成网格
x_min=4
x_max=10
y_min=2.5
y_max=7.5

g1, g2 = np.mgrid[x_min:x_max:200j, y_min:y_max:200j]  # 生成网格采样点
grid_plane = np.stack((g1.flat, g2.flat), axis=1)     # 生成网格，每个点一个坐标
# ----------svm训练开始-----------
# 核：线性，
clf = svm.SVC(kernel='rbf', C=100, gamma=100, probability=True)
# 构造gram的计算方式
#gram = np.dot(X, X.T)
# 训练
clf.fit(X, y)
# ----------svm训练结束-----------


# ----------预测开始-----------
predict=clf.predict(X)
#print(predict)
#print(clf.predict_log_proba(X))
#print(clf.predict_proba(X))
# ----------预测结束-----------

# 几个重要的属性
#print('支撑向量索引：\n',clf.support_)
#print('支撑向量：\n',clf.support_vectors_)
print('每一类支撑向量个数：\n',clf.n_support_)
#print('支撑向量的对偶系数：\n',clf.dual_coef_)
#  print('权重：\n',clf.coef_)     #对非线性，不存在
print('截距：\n',clf.intercept_)
print('0（成功）或者1（失败）：\n',clf.fit_status_)
print('A类概率评估：\n',clf.probA_)
print('B类概率评估：\n',clf.probB_)
print('评估正确率：\n',clf.score(X,y))
# ----------svm训练结束-----------
# 绘制区域网格的
z=clf.predict(grid_plane)
# 控制区域颜色
colors = mpl.colors.ListedColormap([(0.2,0.8,0.8,1), (0.8,0.6,0.5,1)])
# ----------区域网格数据结束-----------

figure=plt.figure(figsize=(8,8))
ax=figure.add_axes([0.1,0.1,0.8,0.8],label='数据集图示')


ax.pcolormesh(g1, g2,z.reshape(g1.shape),cmap=colors)

ax.scatter(x1[0:50],x2[0:50],marker='.',color=(1,0,0,1), label='A类')
ax.scatter(x1[50:100],x2[50:100],marker='x',color=(0,0,1,1), label='B类')

#  绘制支撑向量
ax.scatter(clf.support_vectors_[:clf.n_support_[0],0],clf.support_vectors_[:clf.n_support_[0],1],marker='8',color=(0,1,0,1), label='A类支撑向量', )
ax.scatter(clf.support_vectors_[clf.n_support_[0]:,0],clf.support_vectors_[clf.n_support_[0]:,1],marker='8',color=(1,0,1,1), label='B类支撑向量', )
# 支撑向量平面
# 绘制支撑向量 

ax.set_xbound(x_min,x_max)
ax.set_ybound(y_min,y_max)

plt.legend()
plt.show()

每一类支撑向量个数：
 [46 43]
截距：
 [2.71320865e-05]
0（成功）或者1（失败）：
 0
A类概率评估：
 [-3.56648254]
B类概率评估：
 [0.16228645]
评估正确率：
 0.99

SVM高斯核

惩罚系数C越大，对训练样本集识别率越高，但泛化能力差。惩罚系数越小，训练样本的识别率越低，但相对来说泛化能力强。
一般惩罚系数，按照对数结果方式来取： $10^n$ 。
其中的gamma系数控制分类面的形状， gamma越大，对训练样本集识别率越高，但泛化能力差。惩罚系数越小，训练样本的识别率越低，但相对来说泛化能力强。从上图可以看出gamma越小，越接近线性分类的效果。

逻辑核(sigmoid)

% matplotlib inline
import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
import sklearn as sk
from sklearn import svm
from sklearn import datasets

#加载数据集
data,target = datasets.load_iris(return_X_y=True)
# feature1  feature2  label
X=data[50:150, [0,2]]
x1=X[:,0]
x2=X[:,1]
y=target[0:100]

# 把标签中为0的替换成-1
y[y==0]=-1
# 生成网格
x_min=4
x_max=10
y_min=2.5
y_max=7.5

g1, g2 = np.mgrid[x_min:x_max:200j, y_min:y_max:200j]  # 生成网格采样点
grid_plane = np.stack((g1.flat, g2.flat), axis=1)     # 生成网格，每个点一个坐标
# ----------svm训练开始-----------
# 核：线性，
clf = svm.SVC(kernel='sigmoid', C=100, gamma=0.01, coef0=-1.5,probability=True)
# 构造gram的计算方式
#gram = np.dot(X, X.T)
# 训练
clf.fit(X, y)
# ----------svm训练结束-----------


# ----------预测开始-----------
predict=clf.predict(X)
#print(predict)
#print(clf.predict_log_proba(X))
#print(clf.predict_proba(X))
# ----------预测结束-----------

# 几个重要的属性
#print('支撑向量索引：\n',clf.support_)
#print('支撑向量：\n',clf.support_vectors_)
print('每一类支撑向量个数：\n',clf.n_support_)
#print('支撑向量的对偶系数：\n',clf.dual_coef_)
#  print('权重：\n',clf.coef_)     #对非线性，不存在
#print('截距：\n',clf.intercept_)
print('0（成功）或者1（失败）：\n',clf.fit_status_)
print('A类概率评估：\n',clf.probA_)
print('B类概率评估：\n',clf.probB_)
print('评估正确率：\n',clf.score(X,y))
# ----------svm训练结束-----------
# 绘制区域网格的
z=clf.predict(grid_plane)
# 控制区域颜色
colors = mpl.colors.ListedColormap([(0.2,0.8,0.8,1), (0.8,0.6,0.5,1)])
# ----------区域网格数据结束-----------

figure=plt.figure(figsize=(8,8))
ax=figure.add_axes([0.1,0.1,0.8,0.8],label='数据集图示')


ax.pcolormesh(g1, g2,z.reshape(g1.shape),cmap=colors)

ax.scatter(x1[0:50],x2[0:50],marker='.',color=(1,0,0,1), label='A类')
ax.scatter(x1[50:100],x2[50:100],marker='x',color=(0,0,1,1), label='B类')

#  绘制支撑向量
ax.scatter(clf.support_vectors_[:clf.n_support_[0],0],clf.support_vectors_[:clf.n_support_[0],1],marker='8',color=(0,1,0,1), label='A类支撑向量', )
ax.scatter(clf.support_vectors_[clf.n_support_[0]:,0],clf.support_vectors_[clf.n_support_[0]:,1],marker='8',color=(1,0,1,1), label='B类支撑向量', )
# 支撑向量平面
# 绘制支撑向量 

ax.set_xbound(x_min,x_max)
ax.set_ybound(y_min,y_max)

plt.legend()
plt.show()

每一类支撑向量个数：
 [19 18]
0（成功）或者1（失败）：
 0
A类概率评估：
 [-2.49172675]
B类概率评估：
 [0.13408615]
评估正确率：
 0.94

SVM逻辑（Sigmoid）核

说明

上面的参数的变化规律可以直接调试获取，也可以使用统计表的方式来观察分析，这里也没有做这样的工作。
实际上可以定制SVM核，这里没有参数。

人面猴
序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
沈念sama阅读 214,172评论 6赞 493
死咒
序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
沈念sama阅读 91,346评论 3赞 389
救了他两次的神仙让他今天三更去死
文/潘晓璐我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
开封第一讲书人阅读 159,788评论 0赞 349
道士缉凶录：失踪的卖姜人
文/不坏的土叔我叫张陵，是天一观的道长。经常有香客问我，道长，这世上最难降的妖魔是什么？我笑而不...
开封第一讲书人阅读 57,299评论 1赞 288
港岛之恋（遗憾婚礼）
正文为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
茶点故事阅读 66,409评论 6赞 386
恶毒庶女顶嫁案：这布局不是一般人想出来的
文/花漫我一把揭开白布。她就那样静静地躺着，像睡着了一般。火红的嫁衣衬着肌肤如雪。梳的纹丝不乱的头发上，一...
开封第一讲书人阅读 50,467评论 1赞 292
城市分裂传说
那天，我揣着相机与录音，去河边找鬼。笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。我是一名探鬼主播，决...
沈念sama阅读 39,476评论 3赞 412
双鸳鸯连环套：你想象不到人心有多黑
文/苍兰香墨我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
开封第一讲书人阅读 38,262评论 0赞 269
万荣杀人案实录
序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
沈念sama阅读 44,699评论 1赞 307
护林员之死
正文独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋以下内容为张勋视角年9月15日...
茶点故事阅读 36,994评论 2赞 328
白月光启示录
正文我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
茶点故事阅读 39,167评论 1赞 343
活死人
序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
沈念sama阅读 34,827评论 4赞 337
日本核电站爆炸内幕
正文年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
茶点故事阅读 40,499评论 3赞 322
男人毒药：我在死后第九天来索命
文/蒙蒙一、第九天我趴在偏房一处隐蔽的房顶上张望。院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
开封第一讲书人阅读 31,149评论 0赞 21
一桩弑父案，背后竟有这般阴谋
文/苍兰香墨我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。一阵脚步声响...
开封第一讲书人阅读 32,387评论 1赞 267
情欲美人皮
我被黑心中介骗来泰国打工，没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。一个月前我还...
沈念sama阅读 47,028评论 2赞 365
代替公主和亲
正文我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
茶点故事阅读 44,055评论 2赞 352