Support Vector Machine

支持向量机,一个有监督的机器学习算法。通俗来讲,它是一种二类分类模型,其基本模型定义为特征空间上的间隔最大的线性分类器,其学习策略便是间隔最大化,最终可转化为一个凸二次规划问题的求解。

除了进行线性分类,支持向量机可以使用所谓的核技巧,它们的输入隐含映射成高维特征空间中有效地进行非线性分类。

给定一些数据点,它们分别属于两个不同的类,现在要找到一个线性分类器把这些数据分成两类。如果用 x 表示数据点,用 y 表示类别(y 可以取 1 或者 -1,分别代表两个不同的类),一个线性分类器的学习目标便是要在 n 维的数据空间中找到一个超平面(hyper plane),这个超平面的方程可以表示为
\omega^{T}x+b=0
下面即简单刻画了这些数据:

%matplotlib inline
from sklearn.datasets import make_blobs
import mglearn
import matplotlib.pyplot as plt
#生成聚类数据
X,y = make_blobs(centers=2, random_state=1)
mglearn.discrete_scatter(X[:,0], X[:,1],y)
plt.xlabel("Feature 0")
plt.ylabel("Feature 1")
image.png

根据上图我们可以很容易的作出一条线区分两类数据点.那么要如何确定这个超平面呢?这个超平面应该是 最适合 分开两类数据的直线。而判定“最适合”的标准就是这条直线离直线两边的数据的间隔最大。所以,得寻找有最大间隔的超平面

import numpy as np
from mglearn.plot_helpers import cm2

def plot_2d_separator1(classifier, X, fill=False, ax=None, eps=None, alpha=1,
                      cm=cm2, linewidth=None, threshold=None,
                      linestyle="dashed",a=0):
    # binary?
    if eps is None:
        eps = X.std() / 2.

    if ax is None:
        ax = plt.gca()

    x_min, x_max = X[:, 0].min() - eps, X[:, 0].max() + eps
    y_min, y_max = X[:, 1].min() - eps, X[:, 1].max() + eps
    xx = np.linspace(x_min, x_max, 1000)
    yy = np.linspace(y_min, y_max, 1000)


    X1, X2 = np.meshgrid(xx, yy)

    X_grid = np.c_[X1.ravel(), X2.ravel()-a]

    try:
        decision_values = classifier.decision_function(X_grid)
        levels = [0] if threshold is None else [threshold]
        fill_levels = [decision_values.min()] + levels + [
            decision_values.max()]
    except AttributeError:
        # no decision_function
        decision_values = classifier.predict_proba(X_grid)[:, 1]
        levels = [.5] if threshold is None else [threshold]
        fill_levels = [0] + levels + [1]
    if fill:
        ax.contourf(X1, X2, decision_values.reshape(X1.shape),
                    levels=fill_levels, alpha=alpha, cmap=cm)
    else:
        ax.contour(X1, X2, decision_values.reshape(X1.shape), levels=levels,
                   colors="black", alpha=alpha, linewidths=linewidth,
                   linestyles=linestyle, zorder=5)
from sklearn.svm import LinearSVC
linear_svm = LinearSVC().fit(X,y)
plot_2d_separator1(linear_svm, X,a=2)
plot_2d_separator1(linear_svm, X,a=-2)
mglearn.plots.plot_2d_separator(linear_svm, X)
mglearn.discrete_scatter(X[:,0],X[:,1],y)
plt.xlabel("Feature 0")
plt.ylabel("Feature 1")
image.png

我们注意到 |\omega^{T}x+b| 能够表示点 x 到距离超平面的远近(相对远近),又类标记 y 的符号 (-1, +1) 表示不同类别,故我们可以用
\hat{\gamma}=y(\omega^{T}x+b)=yf(x)
表示间隔,称为函数间隔(functional margin)

但是这样定义间隔,如果成比例的改变 \omegab (如将它们改成 2\omega2b),则函数间隔的值 f(x) 却变成了原来的 2 倍(虽然此时超平面没有改变),所以我们要引入几何间隔(geometrical margin)(经计算得)
\gamma=\frac{\omega^{T}x+b}{||\omega||}=\frac{f(x)}{||\omega||}
为了得到\gamma的绝对值,令 \gamma 乘上对应的类别 y,即可得出几何间隔
\tilde{\gamma}=y\gamma=\frac{\hat{\gamma}}{||\omega||}

对一个数据点进行分类,当超平面离数据点的“间隔”越大,分类的确信度也越大,为了使分类的确信度尽量高,需要让所选择的超平面能够最大化这个“间隔”值。而在这个间隔值上的数据点被称为是支持向量(support vector)的,这就是支持向量机的由来。于是最大间隔分类器(maximum margin classifier) 的目标函数定义为:
max\tilde{\gamma}
同时要满足一些条件,根据间隔的定义,有:
y_{i}(\omega^{T}x_{i}+b)=\hat{\gamma_{i}}\geq \hat{\gamma},\quad i=1,..,n

如果令函数间隔 \hat{\gamma} 等于 1,则有 \tilde{\gamma}=\frac{1}{||\omega||}y_{i}(\omega^{T}x_{i}+b)\geq 1,i=1,...,n,从而上述目标函数转换为
max\frac{1}{||\omega||},\qquad s.t.\ y_{i}(\omega^{T}x_{i}+b)\geq 1,i=1,...,n

由于求 \frac{1}{||\omega||}的最大值相当于求\frac{1}{2}||\omega||^2的最小值,所以上述目标函数等价于
min\frac{1}{2}||\omega||^2,\qquad s.t.\ y_{i}(\omega^{T}x_{i}+b)\geq 1,i=1,...,n

因为现在的目标函数是二次的,约束条件是线性的,所以它是一个凸二次规划问题。此外,由于这个问题的特殊结构,还可以通过拉格朗日对偶性(Lagrange Duality)变换到对偶变量的优化问题,即通过求解与原问题等价的对偶问题得到原始问题的最优解,这就是线性可分条件下支持向量机的对偶算法,这样做的有点在于:一者对偶问题往往更容易求解;二者可以自然的引入核函数。进而推广到非线性分类问题

我们通过给每一个约束条件加上一个拉格朗日乘子(Lagrange multiplier) \alpha,定义拉格朗日函数(通过拉格朗日函数将约束条件融合到目标函数里,从而只用一个函数表达式能够清楚表达问题):
\mathcal{L}(\omega,b,\alpha)=\frac{1}{2}||\omega||^2-\sum_{i=1}^{n}\alpha_{i}(y_{i}(\omega^{T}x_{i}+b)-1)
然后令:
\theta(\omega)=\max_{\alpha_{i}\geq 0}\mathcal{L}(\omega,b,\alpha)
如果没有对 \alpha_{i} 的限制,\theta(\omega) 会等于无穷大

容易验证,当所有约束条件都得到满足时,则有 \theta(\omega)=\frac{1}{2}||\omega||^2 ,亦即最初要最小化的量,所以等价于直接最小化 \theta(\omega)

故目标函数变成:
\min_{\omega,b}\theta(\omega)=\min_{\omega,b}\max_{\alpha_{i}\geq 0}\mathcal{L}(\omega,b,\alpha)=p^*
这里用 p^* 表示这个问题的最优值,且和最初的问题是等价的。如果直接求解,那么一上来便得应对 \omegab 两个参数,而 \alpha_{i} 又是不等式约束,这个求解过程不好做。可以把最小的最大的位置变换一下,变成:
\max_{\alpha_{i}\geq 0}\min_{\omega,b}\mathcal{L}(\omega,b,\alpha)=d^*
交换以后的新问题是原始问题的对偶问题,这个问题的最优值用 d^* 来表示。而且有 d^* \leq p^*。在满足某些条件的情况下,两者等价,这个时候就可以用求解对偶问题来间接地求解原始问题。

上面提到 “d^* \leq p^* 在满足某些条件的情况下,两者等价”,这所谓的“满足某些条件”就是要满足 KKT 条件。经过证明(省略),这里的问题是满足 KKT 条件的,因此现在转化为求解第二个问题。而求解这个对偶问题,分为3个步骤:1.首先要让 \mathcal{L}(\omega,b,\alpha) 关于 \omegab 最小化,2.然后求对 \alpha 的极大,3.最后利用 SMO 算法求解对偶问题中的拉格朗日乘子。

1.首先固定 \alpha ,要让 \mathcal{L} 关于 \omegab 最小化,我们分别对 \omegab 求偏导数,即令 \partial \mathcal{L}/ \partial \omega\partial \mathcal{L}/ \partial b 等于零:
\frac{\partial \mathcal{L}}{\partial \omega}=0 \Rightarrow \omega=\sum_{i=1}^{n}\alpha_{i}y_{i}x_{i}
\frac{\partial \mathcal{L}}{\partial b}=0 \Rightarrow \sum_{i=1}^{n}\alpha_{i}y_{i}=0
将其带入之前的 \mathcal{L},得到:
\mathcal{L}(\omega,b,\alpha)=\sum_{i=1}^{n}\alpha_{i}-\frac{1}{2}\sum_{i,j=1}^{n}\alpha_{i}\alpha_{j}y_{i}y_{j}x_{i}^T x_{j}

2.对 \alpha 求极大,即是关于对偶问题的最优化问题。根据上式有:
\max_{\alpha}\sum_{i=1}^{n}\alpha_{i}-\frac{1}{2}\sum_{i,j=1}^{n}\alpha_{i}\alpha_{j}y_{i}y_{j}x_{i}^T x_{j}
s.t. \quad \alpha_{i}\geq 0,i=1,...,n
\sum_{i=1}^{n}\alpha_{i}y_{i}=0
如果求出了 \alpha_{i} 根据:
\omega^*=\sum_{i=1}^{n}\alpha_{i}y_{i}x_{i}
b^*=-\frac{max_{i:y_{i}=-1}\omega^{*T} x_{i}+min_{i:y_{i}-1}\omega^{*T} x_{i}}{2}
即可求出 \omega,b,最终得出分类超平面和分类决策函数

3.求解上式
\max_{\alpha}\sum_{i=1}^{n}\alpha_{i}-\frac{1}{2}\sum_{i,j=1}^{n}\alpha_{i}\alpha_{j}y_{i}y_{j}x_{i}^T x_{j}
s.t. \quad \alpha_{i}\geq 0,i=1,...,n
\sum_{i=1}^{n}\alpha_{i}y_{i}=0
等价于求解(这个等价引入了后面的松弛变量,及对偶原则,不加证明):
\min_{\alpha}\Psi(\vec \alpha)=\min_{\alpha}\frac{1}{2}\sum_{i=1}^{n}\sum_{j=1}^{n}y_{i}y_{j}K(x_{i},x_{j})\alpha_{i}\alpha_{j}-\sum_{i=1}^{n}\alpha_{i}
s.t. \quad 0\leq \alpha_{i}\leq C,i=1,...,n
\sum_{i=1}^{n}\alpha_{i}y_{i}=0
下面要解决的问题就是:在 \alpha_{i}={\alpha_{1},\alpha_{2},...,\alpha_{n}} 上求上述目标函数的最小值。为了求解这些乘子,每次从中任意抽取两个乘子 \alpha_{1}\alpha_{2},然后固定 \alpha_{1}\alpha_{2} 以外的其它乘子 {\alpha_{3},...,\alpha_{n}},使得目标函数只是关于 \alpha_{1}\alpha_{2} 的函数。这样,不断的从一堆乘子中任意抽取两个求解,不断的迭代求解子问题,最终达到求解原问题的目的。

(具体如何求解,略)

根据前面推导的:
\omega^*=\sum_{i=1}^{n}\alpha_{i}y_{i}x_{i}
带入分类函数有:
f(x)=(\sum_{i=1}^{n}\alpha_{i}y_{i}x_{i})^T x+b=\sum_{i=1}^{n}\alpha_{i}y_{i}<x_{i},x>+b
通过这种形式我们观察到,对于新点x的预测,只需要计算它与训练数据点的内积即可,这一点至关重要,是之后使用Kernel进行非线性推广的基本前提。此外,所谓Supporting Vector 也在这里显示出来——事实上,所有非Supporting Vector 所对应的系数 \alpha 都是等于零的,因此对于新点的内积计算实际上只要针对少量的“支持向量”而不是所有的训练数据即可。

到目前为止,我们的SVM还比较弱,只能处理线性的情况,不过,在得到对偶形式之后,通过Kernel推广到非线性的情况就变成一件比较容易的事情。

我们观察一下下面的数据:

X,y = make_blobs(centers=4, random_state=8)
y = y%2

mglearn.discrete_scatter(X[:,0],X[:,1],y)
plt.xlabel('Feature 0')
plt.ylabel('Feature 1')
image.png

对于上面的数据,我们可以知道他们是线性不可分的,对于非线性的情况,SVM的处理方法是选择一个核函数K(.,.),通过将数据映射到高维空间,来解决在原始空间中先行不可分的问题。

为了理解这个过程(映射到高维空间),我们简单的假设添加第二个特征的平方作为一个新特征,生成一个三维图。

X_new = np.hstack([X, X[:,1:]** 2])
from mpl_toolkits.mplot3d import Axes3D,axes3d
figure = plt.figure()

ax=Axes3D(figure, elev=-152, azim=-26)
mask = y == 0
ax.scatter(X_new[mask,0],X_new[mask,1],X_new[mask,2],c='b',
          cmap=mglearn.cm2, s=60)
ax.scatter(X_new[~mask,0],X_new[~mask,1],X_new[~mask,2],c='r',marker='^',cmap=mglearn.cm2, s=60)
ax.set_xlabel("feature0")
ax.set_ylabel("feature1")
ax.set_zlabel("feature1**2")
image.png

现在我们可以用线性模型将这两个类别分开。

linear_svm_3d = LinearSVC().fit(X_new, y)
coef, intercept = linear_svm_3d.coef_.ravel(),linear_svm_3d.intercept_
figure = plt.figure()
ax = Axes3D(figure, elev=-152, azim=-26)
xx = np.linspace(X_new[:,0].min()-2,X_new[:,0].max()+2,50)
yy = np.linspace(X_new[:,1].min()-2,X_new[:,1].max()+2,50)

XX,YY = np.meshgrid(xx,yy)
ZZ = (coef[0] * XX + coef[1] * YY + intercept) / -coef[2]
ax.plot_surface(XX,YY,ZZ,rstride=8,cstride=8,alpha=0.6)
ax.scatter(X_new[mask,0], X_new[mask,1], X_new[mask,2],c='b',
          cmap=mglearn.cm2,s=60)
ax.scatter(X_new[~mask,0], X_new[~mask,1], X_new[~mask,2],c='r',marker='^',
          cmap=mglearn.cm2,s=60)
image.png

如果将线性SVM模型看作原始特征的函数,那么它实际上已经不是线性的了。他不是一条直线,而是一个椭圆

ZZ = YY**2
dec = linear_svm_3d.decision_function(np.c_[XX.ravel(),YY.ravel(),ZZ.ravel()])
plt.contourf(XX, YY, dec.reshape(XX.shape), levels=[dec.min(),0,dec.max()],cmap=mglearn.cm2, alpha=0.5)
mglearn.discrete_scatter(X[:,0],X[:,1],y)
plt.xlabel("Feature 0")
plt.ylabel("Feature 1")
image.png

从上面的例子中,我们简单的理解了一下处理非线性的方法

支持向量机首先在低维空间中完成计算,然后通过核函数将输入空间映射到高维特征空间,最终在高维特征空间中构造出最优分离超平面,从而把平面上本身不好分的非线性数据分开。

而在我们遇到核函数之前,如果用原始的方法,那么在用线性学习器学习一个非线性关系,需要选择一个非线性特征集,并且将数据写成新的表达形式,这等价于应用一个固定的非线性映射,将数据映射到特征空间,在特征空间中使用线性学习器,因此,考虑的假设集是这种类型的函数:
f(x)=\sum_{i=1}^{N}\omega_{i}\phi_{i}(x)+b
这里的 \phi: \mathcal{X}\to \mathcal{F} 是从输入空间到某个特征空间的映射,这意味着建立非线性学习器分为两步:

1.首先使用一个非线性映射将数据变换到一个特征空间 \mathcal{F};

2.然后在特征空间使用线性学习器分类。

根据之前导出的内积形式:
f(x)=\sum_{i=1}^{n}\alpha_{i}y_{i}<x_{i},x>+b
及上面所述,我们可以在原始特征空间中直接计算内积 <\phi(x_{i}),\phi(x)>,就有可能将两个步骤融合到一起建立一个非线性的学习器,这样直接计算法的方法称为核函数方法

核(Kernel)是一个函数K,对所有的 x,z\in \mathcal{X},满足 K(x,z)=<\phi(x_{i}),\phi(x)>,这里 \phi 是从 \mathcal{X} 到内积特征空间 \mathcal{F} 的映射

an1 = np.linspace(0, 2*np.pi, 100) 
an2 = np.linspace(0, 2*np.pi, 100)
fig,ax = plt.subplots()
ax.plot(3*np.cos(an1)+np.random.rand(100)*0.3, 3*np.sin(an1)+np.random.rand(100)*0.3,'.')
ax.plot(2*np.cos(an2)+np.random.rand(100)*0.3, 2*np.sin(an2)+np.random.rand(100)*0.3,'.')
image.png

对于上面的数据,(用两个椭圆加了点噪声生成),理想的分界应该是一个“圆圈”,我们知道任意一条二次曲线都可以写成下面的形式:
a_{1}X_{1}+a_{2}X_{1}^{2}+a_{3}X_{2}+a_{4}X_{2}^{2}+a_{5}X_{1}X_{2}+a_{6}=0

注意上面的形式,如果我们构造另外一个五维的空间,其中五个坐标的值分别为 Z_{1}=X_{1}, Z_{2}=X_{1}^{2}, Z_{3}=X_{2}, Z_{4}=X_{2}^{2}, Z_{5}=X_{1}X_{2},那么上面的方程在新的坐标系下可以写成:
\sum_{i=1}^{5}a_{i}Z_{i}+a_{6}=0
此时的坐标 Z,是一个超平面方程,也就是,我们做的一个映射 \phi\mathfrak{R}^2 \to \mathfrak{R}^5,将 X 按照上面的规则映射为 Z,那么在新的空间中原来的数据将变成线性可分的(可以参照上面的特征平方的例子)。这正是Kernel方法处理非线性问题的基本思想。

核函数相当于把原来的分类函数:
f(x)=\sum_{i=1}^{n}\alpha_{i}y_{i}<x_{i},x>+b
映射成:
f(x)=\sum_{i=1}^{n}\alpha_{i}y_{i}<\phi(x_{i}),\phi(x)>+b
这样看似我们的问题得到了解决,其实不然,因为我们发现对一个二维空间做映射,得到一个五维空间。如果原始空间是三维的话,那么就会得到一个19维空间,这会给计算带来非常大的困难。如果遇到无穷维的情况,就变得无从计算。

我们从刚才的例子出发。设两个向量 x_{1}=(\eta_{1},\eta_{2})^Tx_{2}=(\xi_{1},\xi_{2})^T 映射后的内积为:
<\phi(x_{i}),\phi(x)>=\eta_{1}\xi_{1}+\eta_{1}^{2}\xi_{1}^{2}+\eta_{2}\xi_{2}+\eta_{2}^{2}\xi_{2}^{2}+\eta_{1}\eta_{2}\xi_{1}\xi_{2}
另外我们注意到:
(<x_{1},x_{2}>+1)^2=2\eta_{1}\xi_{1}+\eta_{1}^{2}\xi_{1}^{2}+2\eta_{2}\xi_{2}+\eta_{2}^{2}\xi_{2}^{2}+2\eta_{1}\eta_{2}\xi_{1}\xi_{2}+1
两者有很多相似的地方,实际上,上面这个式子的计算结果实际上和映射:
\phi(x_{1},x_{2}) =(\sqrt{2}x_{1},x_{1}^{2},\sqrt{2}x_{2},x_{2}^{2},\sqrt{2}x_{1}x_{2},1)^T
之后的内积 <\phi(x_{1}),\phi(x_{2})> 的结果是相等的,它们的区别在于:

1.一个是映射到高维空间中,然后再根据内积的公式进行计算;

2.而另一个则直接在原来的低维空间中进行计算,而不需要显式地写出映射后的结果。

我们把计算两个向量在隐式映射过后的空间中的内积的函数叫做核函数。

常用的核函数有:

1.多项式核 K(x_{1},x_{2})=(<x_{1},x_{2}>+R)^d

2.高斯核 K(x_{1},x_{2})=exp(-||x_{1}-x_{2}||^2 / 2 \sigma^2)(通过调控 \sigma,高斯核实际上具有相当高的灵活性)

3.线性核 K(x_{1},x_{2})=<x_{1},x_{2}>

到此为止,还有些问题没有解决。例如可能并不是因为数据本身是非线性结构的,而只是因为数据有噪声。对于这种偏离正常位置很远的数据点,我们称之为outlier,outlier的存在有可能造成很大的影响,因为超平面本身就是只有少数几个support vector组成的,这样就会造成很大的影响。

为此我们引入松弛变量 \xi_{i},那么约束条件变成:
y_{i}(\omega^T x_{i}+b)\geq 1-\xi_{i},\quad i=1,...,n
它为对应数据点 x_{i} 允许偏移的量,但是 \xi_{i} 又不能随意扩大,所以我们在目标函数后加一项,使得 \xi_{i} 总和最小:
min \frac{1}{2}||\omega||^2 + C\sum_{i=1}^{n}\xi_{i}

其中C是一个参数,用于控制目标函数中两项(“寻找margin最大的超平面”和“保证数据点偏移量最小”)之间的权重。用之前的方法将限制约束条件加入到目标函数中,得到新的拉格朗日函数,如下所示:
\mathcal{L}(\omega,b,\xi,\alpha,r)=\frac{1}{2}||\omega||^2+C\sum_{i=1}{n}\xi_{i}-\sum_{}^{}\alpha_{i}(y_{i}(\omega^T x_{i}+b)-1+\xi_{i})-\sum_{i=1}^{n}r_{i}\xi_{i}
经过计算,现在问题写作:
\max_{\alpha}\sum_{i=1}^{n}\alpha_{i}-\frac{1}{2}\sum_{i,j=1}^{n}\alpha_{i}\alpha_{j}y_{i}y_{j}x_{i}^T x_{j}
s.t. \quad 0 \leq \alpha_{i} \leq C,i=1,...,n
\sum_{i=1}^{n}\alpha_{i}y_{i}=0

这样一个稍微完整的,可以处理线性核非线性并能容忍噪声和outliers的支持向量机差不多简单介绍完了。

©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 199,830评论 5 468
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 83,992评论 2 376
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 146,875评论 0 331
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 53,837评论 1 271
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 62,734评论 5 360
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 48,091评论 1 277
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 37,550评论 3 390
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 36,217评论 0 254
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 40,368评论 1 294
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,298评论 2 317
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,350评论 1 329
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 33,027评论 3 315
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 38,623评论 3 303
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,706评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 30,940评论 1 255
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 42,349评论 2 346
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 41,936评论 2 341

推荐阅读更多精彩内容