九大学习算法简介，玩转分类回归（python版）

机器学习资料多如牛毛，不乏西瓜书，统计学习方法等经典。但多数书籍凸显细节，以致掩盖算法本质，只有参数可调。本文试图避免复杂的公式推导和算法无关的代码，用少量公式和代码，展示算法最基本的组成，以便对机器学习有更加深刻的理解。

线性回归

机器学习定义五花八门，归结起来，就是得到变量之间的期望的依赖关系。
欲拟合变量x,y所满足的线性关系$y=\omega x+b+ϵ$
这里的$y$广义上讲是标签， $x$则是特征，$\omega$是参数，$b$是偏置项，$ϵ$是随机噪声，对于线性模型，我们假定它服从高斯分布。
看$y=f(x)$，$x$是样本空间的向量,$y$是输出空间的标签，我们通常会写成$f(x;\theta )$的形式，这本质就是泛函，参数$\theta$与输入$x$相关，可以看成x的函数，从泛函的角度研究极值，这是优化的基础理论。这给我们一个启示，原始数据至关重要，数据选得好，参数不用调。
解法就是直接求解或是梯度下降法，这里略去过程。
直接得出$\omega =(X^{T}X{)}^{-1}{X}^{T}y$

• 下面展示线性模型的基本代码

#线性模型
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt

import warnings
warnings.filterwarnings(action="ignore", message="^internal gelsd")
mpl.rc('axes',labelsize=14)
mpl.rc('xtick',labelsize=12)
mpl.rc('ytick',labelsize=12)

np.random.seed(24)#固定随机数种子

#直线y = 4x + 3，随机生成100个点
X = 2 * np.random.rand(100,1)
y = 4 + 3 * X + np.random.randn(100,1)
plt.plot(X,y,'b.')
plt.xlabel('$x$',fontsize=18)
plt.ylabel("$y$", rotation=0, fontsize=18)
plt.axis([0, 2, 0, 15])
plt.show()

#正规方程得到斜率和截距，两点确定直线
X_b = np.c_[np.ones((100, 1)), X]  
theta_best = np.linalg.inv(X_b.T.dot(X_b)).dot(X_b.T).dot(y)
X_new = np.array([[0], [2]])
X_new_b = np.c_[np.ones((2, 1)), X_new]  
y_predict = X_new_b.dot(theta_best)

plt.plot(X_new, y_predict, "r-")
plt.plot(X, y, "b.")
plt.axis([0, 2, 0, 15])
plt.xlabel('$x$',fontsize=18)
plt.ylabel("$y$", rotation=0, fontsize=18)
plt.show()
print(theta_best) #得到斜率和截距

#导入机器学习库求解
from sklearn.linear_model import LinearRegression
lin_reg = LinearRegression()
lin_reg.fit(X, y)
print(lin_reg.intercept_, lin_reg.coef_)#结果一样

这段代码的注意事项：随机数种子要固定，不然每次结果都不一样，不太容易对比。我们自定seed时，其实早已经运行过数遍程序，为什么seed不选18，32，因为可能结果更差。

KNN算法

思想非常简单，物以类聚，人以群分。主要难点在如何定义距离，我们不能认为只有欧式空间的距离，实际上距离的定义就有很多种。距离的特性，非负性，对称性，三角不等式，只要满足这三个条件即可。k 值的选择反映了对近似误差与估计误差之间的权衡，通常由交叉验证选择最优的 k 。

• 代码如下

from sklearn import datasets
from sklearn.model_selection import train_test_split,cross_val_score
import matplotlib.pyplot as plt
from sklearn.neighbors import KNeighborsClassifier

#载入数据集
iris = datasets.load_iris()
X = iris.data
y = iris.target

train_X,test_X,train_y,test_y = train_test_split(X,y,test_size=0.2,random_state=3)

#搜索最佳K值
cv_scores =[]
for n in range(1,31):
    knn = KNeighborsClassifier(n)
    scores = cross_val_score(knn,train_X,train_y,cv=10,scoring='accuracy')
    cv_scores.append(scores.mean())

plt.plot(range(1,31),cv_scores)
plt.xlabel('K')
plt.ylabel('scores')
plt.show()

这段代码的注意事项：train_test_split,cross_val_score这两个函数有划分数据集和评估模型的作用，详情请查看sklearn文档。

支持向量机

利用了最大超平面间隔的思想，本质上是二次规划问题，即
$$min\frac{1}{2}\parallel \omega {\parallel }^2$$$$s.t.y_i(w^T{x}_i+b)\ge 1,i=1,2,...,m$$用拉格朗日乘子法和KKT条件求解最优值，转化成对偶问题的形式，过程略。
对偶问题
$$\min \frac{1}{2}\sum _{i=1}^N\sum _{j=1}^N{\alpha }_i{\alpha }_j{y}_i{y}_j\left(x_i\cdot x_j\right)-\sum _{i=1}^N{\alpha }_i$$

$$s.t.\sum _{i=1}^N{\alpha }_i{y}_i=0$$

$${\alpha }_i⩾0,i=1,2,\cdots ,N$$

可采用SMO算法求解，特点是不断地将原二次规划问题分解为只有两个变量的二次规划子问题，并对子问题进行解析求解，直到所有变量满足KKT条件为止。

支持向量可在间隔边界上，也可在间隔边界与分离超平面之间，或者在分离超平面误分一侧。最优分离超平面由支持向量完全决定。

线性支持向量机学习等价于最小化二阶范数正则化的合页函数

$$\sum _{i=1}^N{\left[1-y_i\left(w\cdot x_i+b\right)\right]}_{+}+\lambda \Vert w{\Vert }^2$$
对于非线性分类，引入核函数，$K(x,z)$是一个核函数，或正定核，意味着存在一个从输入空间x到特征空间的映射$\mathcal{X}\to \mathcal{H}$，对任意$\mathcal{X}$，有

$$K(x,z)=\varphi (x)\cdot \varphi (z)$$

• 代码如下

import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasets
from sklearn.svm import LinearSVC
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import PolynomialFeatures, StandardScaler
from matplotlib.colors import ListedColormap

#产生数据点
X, y = datasets.make_moons(noise=0.15, random_state=666)

def PolynomialSVC(degree, C=1.0):
    return Pipeline([
        ("poly", PolynomialFeatures(degree=degree)),
        ("std_scaler", StandardScaler()),
        ("linearSVC", LinearSVC(C=C))
    ])


poly_svc = PolynomialSVC(degree=3)
poly_svc.fit(X, y)

#绘制决策边界
def plot_decision_boundary(model, axis):

    x0, x1 = np.meshgrid(
        np.linspace(axis[0], axis[1], int((axis[1] - axis[0])*100)).reshape(1, -1),
        np.linspace(axis[2], axis[3], int((axis[3] - axis[2])*100)).reshape(1, -1)
    )
    X_new = np.c_[x0.ravel(), x1.ravel()]

    y_predict = model.predict(X_new)
    zz = y_predict.reshape(x0.shape)

    custom_cmap = ListedColormap(['#EF9A9A', '#FFF59D', '#90CAF9'])
    plt.contourf(x0, x1, zz, cmap=custom_cmap)

plot_decision_boundary(poly_svc, axis=[-1.5, 2.5, -1.0, 1.5])
plt.scatter(X[y == 0, 0], X[y == 0, 1])
plt.scatter(X[y == 1, 0], X[y == 1, 1])
plt.xlabel('x')
plt.ylabel('y')
plt.show()

注意事项：make_moons双月数据集的参数，Pipeline管道流的用法，起到并行加速的作用。核函数的选取，有多项式核，高斯核，拉普拉斯核，sigmoid核等。

决策树

决策树可以转换成一个if-then规则的集合，也可以看作是定义在特征空间划分上的类的条件概率分布。决策树由多个阶跃函数或符号函数组合而成，可以逼近十分复杂的函数。
算法有ID3，C4.5，CART等，基于信息论原理，重点把握信息增益的概念。
终止条件有如下几种：
1、特征已经用完了：没有可供使用的特征再进行分裂了，则树停止分裂；
2、子节点中没有样本了：此时该结点已经没有样本可供划分，该结点停止分裂；
3、树达到了人为预先设定的最大深度：depth >= max_depth，树停止分裂。
4、节点的样本数量达到了人为设定的阈值：样本数量 < min_samples_split ，则树停止分裂。
决策树剪枝方法防止过拟合，也需要关注。

Bagging

方法如下,设大小为 $N$ 的样本$X$。我们可以通过随机均匀的从训练集中抽取 $N$ 个元素，作为采集池，有放回重复进行N次。假设，我们有一个训练集 $X$，利用自举，我们生成样本 $X_1,X_2,\dots ,X_M$ 。现在，为每一个自举的样本，我们都去训练一个专属于它的分类器 $a_i(x)$ 。最终的分类器是平均所有这些单独分类器的输出。在分类算法中，这种技术称之为投票。

• 代码如下

from sklearn.tree import DecisionTreeClassifier
from sklearn.metrics import accuracy_score
from sklearn.model_selection import train_test_split
from sklearn.datasets import make_moons
import matplotlib.pyplot as plt
import numpy as np

X, y = make_moons(n_samples=500, noise=0.2, random_state=24)
X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=8)

tree_clf = DecisionTreeClassifier(random_state=24)
tree_clf.fit(X_train, y_train)
y_pred_tree = tree_clf.predict(X_test)
print(accuracy_score(y_test, y_pred_tree))

from matplotlib.colors import ListedColormap

def plot_decision_boundary(clf, X, y, axes=[-1.5, 2.5, -1, 1.5], alpha=0.5, contour=True):
    x1s = np.linspace(axes[0], axes[1], 100)
    x2s = np.linspace(axes[2], axes[3], 100)
    x1, x2 = np.meshgrid(x1s, x2s)
    X_new = np.c_[x1.ravel(), x2.ravel()]
    y_pred = clf.predict(X_new).reshape(x1.shape)
    custom_cmap = ListedColormap(['#fafab0','#9898ff','#a0faa0'])
    plt.contourf(x1, x2, y_pred, alpha=0.3, cmap=custom_cmap)
    if contour:
        custom_cmap2 = ListedColormap(['#7d7d58','#4c4c7f','#507d50'])
        plt.contour(x1, x2, y_pred, cmap=custom_cmap2, alpha=0.8)
    plt.plot(X[:, 0][y==0], X[:, 1][y==0], "yo", alpha=alpha)
    plt.plot(X[:, 0][y==1], X[:, 1][y==1], "bs", alpha=alpha)
    plt.axis(axes)
    plt.xlabel(r"$x_1$", fontsize=18)
    plt.ylabel(r"$x_2$", fontsize=18, rotation=0)

from sklearn.ensemble import BaggingClassifier
bag_clf = BaggingClassifier(
    DecisionTreeClassifier(random_state=24), n_estimators=500,
    max_samples=100, bootstrap=True, n_jobs=-1, random_state=24)
bag_clf.fit(X_train, y_train)

plt.figure(figsize=(12,6))
plt.subplot(121)
plot_decision_boundary(tree_clf, X, y)
plt.title("Decision Tree", fontsize=14)
plt.subplot(122)
plot_decision_boundary(bag_clf, X, y)
plt.title("Decision Trees with Bagging", fontsize=14)
plt.show() 

注意事项：这里的BaggingClassifier用决策树来投票。

随机森林

由多颗决策树组成，同理。

• 代码如下

#随机森林
import numpy as np
import matplotlib.pyplot as plt
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import train_test_split
from sklearn.datasets import make_moons
from matplotlib.colors import ListedColormap

np.random.seed(24)#固定随机数种子

X, y = make_moons(n_samples=500, noise=0.1, random_state=24)
X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=24)

def plot_decision_boundary(clf, X, y, axes=[-1.5, 2.5, -1, 1.5], alpha=0.5, contour=True):
    x1s = np.linspace(axes[0], axes[1], 100)
    x2s = np.linspace(axes[2], axes[3], 100)
    x1, x2 = np.meshgrid(x1s, x2s)
    X_new = np.c_[x1.ravel(), x2.ravel()]
    y_pred = clf.predict(X_new).reshape(x1.shape)
    custom_cmap = ListedColormap(['#fafab0','#9898ff','#a0faa0'])
    plt.contourf(x1, x2, y_pred, alpha=0.3, cmap=custom_cmap)
    if contour:
        custom_cmap2 = ListedColormap(['#7d7d58','#4c4c7f','#507d50'])
        plt.contour(x1, x2, y_pred, cmap=custom_cmap2, alpha=0.8)
    plt.plot(X[:, 0][y==0], X[:, 1][y==0], "yo", alpha=alpha)
    plt.plot(X[:, 0][y==1], X[:, 1][y==1], "bs", alpha=alpha)
    plt.axis(axes)
    plt.xlabel(r"$x_1$", fontsize=18)
    plt.ylabel(r"$x_2$", fontsize=18, rotation=0)

rnd_clf = RandomForestClassifier(n_estimators=500, n_jobs=-1, random_state=24)
rnd_clf.fit(X_train,y_train)
plt.figure(figsize=(6, 4))
plot_decision_boundary(rnd_clf,X_test,y_test)

plt.show()

朴素贝叶斯

朴素贝叶斯法是典型的生成学习方法。生成方法由训练数据学习联合概率分布
$P(X,Y)$，然后求得后验概率分布$P(Y|X)$。具体来说，利用训练数据学习$P(X|Y)$和$P(Y)$的估计，得到联合概率分布：

$$P(X,Y)＝P(Y)P(X|Y)$$
将输入$x$分到后验概率最大的类$y$。

$$y=\arg \underset{c_k}{\max }P\left(Y=c_k\right)\prod _{j=1}^{n}P\left(X_j=x^{(j)}|Y=c_k\right)$$

后验概率最大等价于0-1损失函数时的期望风险最小化。

• 代码如下

from sklearn.datasets import make_blobs
from sklearn.model_selection import train_test_split
from sklearn.naive_bayes import GaussianNB
import matplotlib.pyplot as plt
import numpy as np
np.random.seed(24)
# 产生数据集
X, y = make_blobs(n_samples=200, n_features=3, centers=4, random_state=8)  # X是数据特征值集合，y是类别
X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=8)
#根据X来划分成两类
clf = GaussianNB()
clf.fit(X_train, y_train)
plt.scatter(X[X[:,0]<4][:,0],X[X[:,0]<4][:,1])
plt.scatter(X[X[:,0]>4][:,0],X[X[:,0]>4][:,1])
plt.xlabel('X')
plt.ylabel('y')
#绘图直接选取两点
y_pred=clf.predict(X_test)
score=clf.score(X_test, y_test)
print(score)
plt.plot([2,6],[10,-7.5],'r-')
plt.xlabel('X')
plt.ylabel('y')
plt.show()

注意事项：只需要牢牢把握基本代码框架即可，对于多分类问题，则需要转化为二分类问题。

K均值聚类

聚类是针对给定的样本，依据它们属性的相似度或距离，将其归并到若干个“类”或“簇”的数据分析问题。一个类是样本的一个子集。直观上，相似的样本聚集在同类，不相似的样本分散在不同类。步骤如下：

第一步，随机选择 K 个点, 或者指定
第二步，遍历所有点 P，将 P 放入最近的聚类中心的集合中
第三步，遍历每一个点集，计算新的聚类中心，重复直到收敛

• 层次聚类

聚合（自下而上）：聚合法开始将每个样本各自分裂到一个类，之后将相距最近的两类合并，建立一个新的类，重复次操作知道满足停止条件，得到层次化的类别。

分裂（自上而下）： 分裂法开始将所有样本分到一个类，之后将已有类中相距最远的样本分到两个新的类，重复此操作直到满足停止条件，得到层次化的类别。

• 代码如下

from sklearn.datasets import make_blobs 
import matplotlib.pyplot as plt
import numpy as np
np.random.seed(24)
# 产生数据集
X, y = make_blobs(n_samples=1000, n_features=2, centers=4, random_state=8)  # X是数据特征值集合，y是类别
 
fig, ax1 = plt.subplots(1)  # 创建一个子图,返回 Figure对象（fig) 和 子图对象（ax1）
ax1.scatter(X[:, 0], X[:, 1], marker='o', s=8) 
plt.show()

# 不同簇显示不同的颜色，这里是数据实际的分类
color = ["red","pink","orange","gray"] 
fig, ax1 = plt.subplots(1)
for i in range(4):
    ax1.scatter(X[y==i, 0], X[y==i, 1], marker='o', s=8, c=color[i])
plt.xlabel('X')
plt.ylabel('y')
plt.show()

from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_score
from sklearn.metrics import silhouette_samples 
K = [3,4,5,6]  # 猜测簇的个数
# KMeans对象有三个属性：labels_,cluster_centers_,inertia_
# 属性 labels_：模型 聚类 得到的类别（每个样本所对应类）
# 属性 cluster_centers_：最终得到的所有的质心
# 属性 inertia_：总距离平方和---越小越好
scores = []
for k in K:
    cl = KMeans(n_clusters=k,random_state=24).fit(X)
    print(cl.inertia_)
    score = silhouette_score(X, cl.labels_)
    scores.append(score)
#实际上，当我们设的 k 越大（不超过样本数），总距离平方和就会越小。这是因为，k越大，则质心就越多，当每一个数据点都是一个质心，那此时总距离平方和就会等于0

print()
print(scores)

注意事项：我们用总距离平方和，轮廓系数来评价聚类效果。当我们所设的簇越大（不超过样本数），总距离平方和就会越小。轮廓系数会结合内聚度和分离度两种因素来考虑。

主成分分析

需要结合奇异值分解和因子分析来看，要求复习矩阵相关知识。奇异值是特征值的平方根，任何矩阵都可以奇异值分解，实对称矩阵可以特征分解，因子分析通过特征组合创造新的因子，具有实际意义。

假设$x$为$m$ 维随机变量，其均值为$\mu$，协方差矩阵为$\Sigma$。

考虑由$m$维随机变量$x$到$m$维随机变量$y$的线性变换
$$y_i={\alpha }_i^{T}x=\sum _{k=1}^m{\alpha }_{ki}{x}_k,i=1,2,\cdots ,m$$

其中${\alpha }_i^T=({\alpha }_{1i},{\alpha }_{2i},\cdots ,{\alpha }_{mi})$。

如果该线性变换满足以下条件，则称之为总体主成分：

（1）${\alpha }_i^T{\alpha }_i=1,i=1,2,\cdots ,m$；

（2）$\mathrm{cov}(y_i,y_j)=0(i\ne j)$;

（3）变量$y_1$是$x$的所有线性变换中方差最大的；$y_2$是与$y_1$不相关的$x$的所有线性变换中方差最大的；一般地，$y_i$是与$y_1,y_2,\cdots ,y_{i-1},(i=1,2,\cdots ,m)$都不相关的$x$的所有线性变换中方差最大的；这时分别称$y_1,y_2,\cdots ,y_m$为$x$的第一主成分、第二主成分、…、第$m$主成分。
矩阵$X$的奇异值分解算法。针对$m\times n$样本矩阵$X$
$$X^{\prime }=\frac{1}{\sqrt{n-1}}{X}^T$$对矩阵$X^{\prime }$进行截断奇异值分解，保留$k$个奇异值、奇异向量，得到
$$X^{\prime }=US{V}^T$$ $V$的每一列对应一个主成分，得到$k\times n$样本主成分矩阵$Y$ $$Y=V^{T}X$$

• 代码如下

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.datasets import load_iris
from sklearn.preprocessing import StandardScaler

np.random.seed(24)
#加载数据
iris = load_iris()
X = iris.data
X_norm = StandardScaler().fit_transform(X)
X_norm.mean(axis=0)

# 求特征值和特征向量
ew,ev = np.linalg.eig(np.cov(X_norm.T))
# 特征向量特征值的排序
ew_oreder = np.argsort(ew)[::-1]
ew_sort = ew[ew_oreder]
ev_sort = ev[:, ew_oreder]  # ev的每一列代表一个特征向量
ev_sort.shape # (4,4)

# 我们指定降成2维， 然后取出排序后的特征向量的前两列就是基
K = 2
V = ev_sort[:, :2]  # 4*2

# 最后，我们得到降维后的数据
X_new = X_norm.dot(V)    # shape (150,2)

colors = ['red', 'black', 'orange']

plt.figure()
for i in [0, 1, 2]:
    plt.scatter(X_new[iris.target==i, 0],
                X_new[iris.target==i, 1],
                alpha=.7,
                c=colors[i],
                label=iris.target_names[i]
               )

plt.legend()
plt.title('PCa of IRIS dataset')
plt.xlabel('PC_0')
plt.ylabel('PC_1')
plt.show()

#使用库
from sklearn.decomposition import PCA

# 然后使用
pca = PCA(n_components=2)
X_new = pca.fit_transform(X_norm)

"""查看PCA的一些属性"""
print(pca.explained_variance_)    # 属性可以查看降维后的每个特征向量上所带的信息量大小（可解释性方差的大小）
print(pca.explained_variance_ratio_)  # 查看降维后的每个新特征的信息量占原始数据总信息量的百分比
print(pca.explained_variance_ratio_.sum())    # 降维后信息保留量
pca_line = PCA().fit(X_norm)
plt.plot([1,2,3,4],np.cumsum(pca_line.explained_variance_ratio_))
plt.xlabel('component')
plt.ylabel('proportion')
plt.xticks([1,2,3,4])
plt.show()

注意事项：对于Iris数据集，前2个主成分已经大于0.95，说明这个数据集太容易区分，这也是为什么书上都会用这个例子，换别的数据集可能效果很差，不具有说服力，还以为是算法的问题。

总结

相比厚厚的一本书，本文可谓极其简洁，忽略了很多知识点，理解起来可能比较困难，基本回顾了机器学习的主要方法，然而举例比较简单，主要是针对代码不太熟练的同学而写，可以通过运行一个简单的示例来加深理解，实际上，可视化是非常困难的，所举例子都是二维情况，显得直观，而把真正的困难隐藏起来，这也是大多书籍常用手段。很多入门书籍对简单的细节重复无数次，而缺乏深入思考，不过是文档教程的翻译和搬运，这是需要极力避免的。希望本文能对大家学习机器学习有所帮助。
代码地址