模型评估与选择

2.1经验误差与过拟合

经验误差(empirical error)为训练误差(training error)，即在训练集上得到的误差，
在新样本上的误差为泛化误差(generalization error)

当学习器把训练样本学习的太好时，会产生过拟合(overfitting)，于此相对的是欠拟合(underfitting)即对训练样本的一般性质上尚未学好

过拟合最常见的情况是学习能力过于强大，以至于把训练样本所有的不太一般的性质也学会了，而欠拟合通常是由于学习能力低下导致的

因为数据的不确定和不准确性，所以没法得到一个完全正确的模型，所以过拟合不可避免。即因为数据的随机性和其他原因，经验误差等于0，但是泛化误差也等于0的情况基本上是不存在的

2.2 评估方法

通过实验来得到的泛化误差来评估并进行选择
通常需要一个测试集(testing set)来测试学习器对新样本的判断能力，用测试误差(testing error)作为泛化误差的一个近似

对于一个m个样本的数据集D，得到一份训练集S和一份测试集T

2.2.1 留出法

即将数据集分成两份数据集，当做测试集和训练集，这两个数据集的样本都不一样

需要注意的是S和T的分布需要一致，避免划分是引入额外的误差，列如分类中，各分类的比例保持一致

从采样(sampling)的角度就是分层采样(stratified sampling)，即按照类别分别采集相同比例的样本

假设D有三个类别，分别占有4，4，2的比例，那假设训练样本抽取80%的数据，则对于占40%的class1取其中的80%，同理对其他的类别

但是留出法得到的结果往往不稳定，因为根据采样的方法不同，得到的D,T样本也不同，因为如果对D进行排序，则取前面的80%数据和取后面80%的数据，得到的不一样，
一般采用若干次随机划分、重复进行实验评估后取平均值

但是实际是评估数据集D训练出来的学习器性能，所以常用的比例大约是2/3到4/5的数据用于当做训练集，来保证能尽可能接近于训练集D的分类器，但是又有一定数据来评估，得到较准确的评估值

2.2.2 交叉验证法

将数据集D划分为k个大小类似的互斥子集，每个子集都能保持一样的数据分布，也就是用了分层采样得到。

将k-1个子集当做训练集，1个当做测试集，那就能得到k个评估，取平均值作为整体的评估结果。
常用的k值为10，也就是k-折交叉验证，其他常用的k值为5，20。

其中如果将k定为m,则测试集只是其中一个样本，叫做留一法，这样可以保证得到的学习器可以和D训练出来的近似一致，
但是没法保证得到的估计结果会比其他评估方法好，并且计算量太大

2.2.3 自助法

上面的方法因为都划分了一部分当成测试集，导致实际学习器训练的数据比实际数据D小，所以必然引入了一些因训练样本规模不一样导致的估计偏差，留一法虽然影响较小，但是计算量提高了。

为了解决训练规模的影响，可使用自助法

自助法(bosststrapping)是以自助采样法(boostrap sampling)为基础，给定m个数据集D，对D进行采样产生数据集D’,然后在将该样本放回D中，使得样本在下次采样中可以被采集到,就是放回的采样；重复m次后，得到了m个样本的数据D’；就是对数据D进行放回采样，每次采集m个样本。

这样样本D’的数据规模和原本的一样，并且因为有些数据被重复抽到导致一部分的数据并没有包括

数据在m次采样中始终不包括的概率是$(1-\frac{1}{m})^m$,取极限得到

即通过自助法，还有大概36.8%的样本没有被抽到，可以用D’做训练集，D\D’当测试集，这样的测试结果为包外估计(out-of-bag estimate)
自助法在数据集较小、难以有效划分训练集/测试集时很有效；此外，自助法能从初试数据集中产生不同的训练集，对集成学习等方法有很大的好处，但是因为改变了原本数据的分布，所以引入了估计偏差。因此在初试数据量足够时，留出法和交叉验证法较为常用

2.2.4 调参和最终模型

大部分的算法都有参数(parameter)需要设置，除了对模型算法进行选择，还需要对算法参数进行设置，叫做参数调节即调参(parameter tuning)
许多参数都是在一个范围内选择的，参数往往不是最佳的，这是在计算开销和性能估计做的折中选择。
通常得到的最终模型还需要对整个训练进行训练，用的就是得到的最佳参数和算法模型，因为在之前得到的模型和算法参数都是根据部分数据集训练得到的。
为了和最终测试的测试集加以区分，训练数据中来估计的数据集称之为验证集(validation set)

2.3 性能度量

除了实验估计方法还需要衡量模型泛化能力的评价标准，就是性能度量(performace measure)。在任务需求下，使用不同的性能度量，可以得到不同的结果。所以模型的好坏除了算法和数据，还取决于任务需求。
回归任务最常用的性能度量是 均方误差(mean squared error)

更一般的，对于数据分布$\mathcal{D}$和概率密度函数$\mathcal{p}(\cdot)$，均方误差可描述为

2.3.1 错误率和精度

对于数据集D，分类错误率定义为

精度则为

更一般的，对于数据分布$\mathcal{D}$和概率密度$\mathcal{p}(\cdot)$来说，错误率和精度为

精度则为

2.3.2 查准率、查全率和F1

如果不为了得到准确率，而是关心判断正确的样本中正确的比例，或者 有多少正例的样本被判断对了，则需要其他性能度量方法
查准率(precision)和查全率(recall)亦称为召回率更适合用于这类需求的性能度量

对于二分问题，根据真实类别和学习器判断的类别组合成四个分类

1. 真正例(true positive) 真实为 正，预测也为 正 TP
2. 假正例(false positive) 真实为 负，预测为 正 FP
3. 真反例(true negative) 真实为 负， 预测为 负 TN
4. 假反例(false negative) 真实为 正，预测为 负 FN

查准率P和查全率R分别为

查准率和查全率互相制约，基本上一个比较高时，另一个则较低，因为如果想要查全率高，则需要尽量的判断为正，但是会增加假正例，所以降低了查准率。如果要查准率高，则尽量选择有把握的正例，则导致了真正例在所有正例的比例降低，所以查全率会降低。

P-R曲线
根据学习器对样本的为正例的判断可信度排序，得到从最有可能是正例到最有可能是负例排序的样本，根据每个样本为阀值进行划分，判断前面的为正例，后面的为负例，得到查准率和查全率，用查全率作为横坐标，查准率作为纵坐标，得到曲线就是P-R曲线

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
%matplotlib inline
y = np.array([0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1]) #预测的类别
y_hat = np.array([0,0,0,0,0,1,0,0,1,0,0,1,1,1,1,1,0,1,1,0,1,1,1]) #实际的类别

def f(y,y_hat):
    #默认已经排序了，不写排序的代码
    m = y_hat.shape[0]
    x = []
    y = []
    for i in range(1,m + 1):
        positive = y_hat[-i:]
        nagetive = y_hat[:-i]
        tp = positive.sum()
        fp = len(positive) - tp
        fn = nagetive.sum()
        tn = len(nagetive) - fn
        p = tp/float(len(positive))
        r = tp / float(y_hat.sum())
        x.append(r)
        y.append(p)
    return x,y

res = f(y,y_hat)

plt.plot(res[0],res[1])

[<matplotlib.lines.Line2D at 0x1a1b31d198>]

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

from sklearn.metrics import precision_recall_curve
from sklearn.utils.fixes import signature
plt.figure("P-R Curve")
plt.title('Precision/Recall Curve')
plt.xlabel('Recall')
plt.ylabel('Precision')
#y_true为样本实际的类别，y_scores为样本为正例的概率
y_true = np.array([1, 1, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 0, 0, 0, 0, 1, 0, 0, 0])
y_scores = np.array([0.9, 0.75, 0.86, 0.47, 0.55, 0.56, 0.74, 0.62, 0.5, 0.86, 0.8, 0.47, 0.44, 0.67, 0.43, 0.4, 0.52, 0.4, 0.35, 0.1])
precision, recall, thresholds = precision_recall_curve(y_true, y_scores)
#print(precision)
#print(recall)
#print(thresholds)
plt.plot(recall,precision)
#precision, recall, thresholds = precision_recall_curve(y_hat, y)
#print(recall)
#print(precision)
#plt.plot(recall,precision)

[<matplotlib.lines.Line2D at 0x1a1b2a14a8>]

如果学习器A的P-R曲线包裹学习器B的曲线，则表明A优于B，但是曲线相交，则用曲线下的面积来判断，但是这个值很难判断，所以有其他判断方式

平衡点(Break-Even Point 简称BEP)，即查准率=查全率的时候.
BEP过于简单，所以一般使用F1度量,

通过对查准率和查全率的偏向，得到一个更一般的$F_\beta$,定义为

其中$\beta>0$,度量了查全率对查准率的相对重要程度，当$\beta=1$就是标准的F1,$\beta >1$则说明查全率更重要，$\beta<1$说明查准率更重要

假设进行了多次的训练/测试，就获得了多个二分类的混淆矩阵，要综合判定，则有两种方法

一种是分别计算查准率和查全率，记为(P1,R1),(P2,R2),…$(P_n,R_n)$,在计算平均值，得到宏查准率(macro-P)，宏查全率(macro-R)，得到相应的宏F1(macro-F1)

另一种为对混淆矩阵的元素，求平均值，得到对应的$\overline{TP}$,$\overline{TN}$,$\overline{FP}$,$\overline{FN}$
求得微查准率(micro-P)，微查全率(micro-R)，和对应的微F1(micro-F1)

$$micro-R = \frac{\overline{Tp}}{\overline{TP} +\overline{FN}}$$
$$micro-F1 = \frac{2 \times  micro-P \times  micro-R}{micro-P + micro-R}$$

2.3.3 ROC和AUC

P-R曲线对于样本不均衡的情况下，不稳定，所以可以使用ROC来判断

ROC(Receiver Operation Characteristic)全称是受试者工作特征曲线，与P-R曲线类似，但是横纵坐标分别是 假正例率(False Positive Rate FPR)和 真正例率(True Positive Rate TPR)分别表示为

如果学习器A的ROC曲线包裹学习器B的曲线，则表明A优于B,但是如果曲线交叉，则用去线下的面积AUC(Area Under ROC Curve)来判断好坏。

假设ROC曲线是由{$(x_1,y_1),(x_2,y_2),…,(x_n,y_n$}的点组成的，则AUC可以估计为

可以看到AUC和样本预测的排序质量有关，假设有$m^+$,$m^-$个正反例，$D^+,D^-$表示正反例集合，则排序的损失(loss)函数为

即考虑每对正反例，若正例预测值小于反例，则记一个”罚分”，若相等，则记半个”罚分”，因为学习器是预测正例的，所以理想的状态是，正例的最小预测值大于反例的最大预测值，$\ell_{rank}$对应的就是ROC曲线之上的面积；

若一个正例在ROC曲线上对应标记点坐标为(x,y),则x恰好是排序在其之前的反例所占的比例，即假正例率，因此有

ROC曲线与P-R曲线的区别

最主要的区别是ROC可以在样本不均衡的情况下做判断，而P-R曲线会有很大的影响，但是P-R曲线更加直观

2.3.4代价敏感错误率与代价曲线

不同哦该的错误造成的后果也不同，为不同类型错误所造成的不同损失，可为错误赋予非均等的代价(unqeual cost).

根据领域知识得到一个代价矩阵(cost matrix),其中$cost_{ij}$表示第i类被预测为j类，一般来说$cost_{ii}$为0

在非均等代价下，我们需要的不再是简单的最小化错误次数，而是最小化总体代价(total cost),将第0类作为正例，第1类作为反例，则代价敏感(cost-sensitive)的错误率为

在非均等代价下，ROC不能直接反应处学习器的期望总体代价，所以需要用代价曲线(cost curve)来达到目的,
横坐标为取值在[0,1]的正例概率代价

其中p为样例为正例的概率，纵坐标为取值为[0,1]的归一化代价

其中FPR为假正率，而FNR=1-TPR为假反例率

2.4 比较检验

统计假设检验(hypothesis test)为我们在学习器性能比较提供了重要的依据，默认以错误率作为性能度量，用$\epsilon$表示

2.4.1 假设检验

假设检验中的”假设”是对学习器泛化错误分布的某种判断或猜测，例”$\epsilon=\epsilon_0$”。在实际任务中，并不知道泛化错误率，只能知道其测试错误率$\hat{\epsilon}$,泛化错误率和测试错误率未必相同，但是两者接近的可能性较大，相差很远的可能性很小。因此可根据测试错误率推出泛化错误率的分布

假设泛化错误为$\epsilon$,测试错误率为$\hat{\epsilon}$,$\hat{\epsilon}$意味着如果有m个样本，则会有$\hat{\epsilon}\times m$个样本会被分错。
假设测试样本，有$m’$个样本被分错,则泛化错误率为$\epsilon$的学习器将$m’$个样本分错的概率为$\tbinom{m}{m’}\epsilon^{m’}(1-\epsilon)^{m-m’}$,因此可以估计出恰好将$\hat{\epsilon}\times m$个样本分错的概率为

对此求导，解$\frac{\partial P(\hat{\epsilon};\epsilon)}{\partial \epsilon}=0$得到$p(\hat{\epsilon};\epsilon)$在$\epsilon=\hat{\epsilon}$时最大，这符合二项分布(binomial)

使用二项检验(binomial test)来对『$\epsilon \leq \epsilon_0$』这样的假设进行验证

在$1-\alpha$的概率内所能得到的最大错误率为

其中$1-\alpha$反映了置信度(confidence),因为二项分布，随着p的变大，会向右偏，所以随着p的变大，会不满足条件，所以有个极值,就是需要得到这个最大值错误率，如果测试错误率小于这个最大错误率，则说明，在$1-\alpha$的置信度里，假设成立的，不能被拒绝，即认为在$1-\alpha$的置信度下，测试错误率为$\hat{\epsilon}$的学习器的泛化错误率不大于$\epsilon_0$。反之则需要拒绝

但是一般不止做一次留出法估计，或者其交叉验证中k取2的估计，会进行很多次训练/测试,所以会得到很多错误率，这时可以用t验证(t-test)。

假设得到了k个测试错误率$\hat{\epsilon_1},\hat{\epsilon_2},…,\hat{\epsilon_k}$,可以求的错误率$\mu$,方法$\sigma^2$为

考虑到k个测试错误率可以看做繁华错误率$\epsilon_0$的独立采样，则t统计量为

服从自由度为k-1的t分布

得到的$\tau_t$如果小于查表得到的t值，则接受原假设，否则拒绝原假设

2.4.2 交叉验证t检验

对于两个学习器A和B，使用k折交叉验证法得到的测试错误率分别为$\epsilon_1^A,\epsilon_2^A,…,\epsilon_k^A$和$\epsilon_1^B,\epsilon_2^B,…,\epsilon_k^B$,其中$\epsilon_i^A$和$\epsilon_i^B$是相同的第i折训练/测试集上得到的结果，则可用k折交叉验证”成对t验证”(paired t-tests)来进行比较。这里的思想为如果两个学习器的性能相同，则它们使用相同的训练/测试集得到的测试错误率应相同，即$\epsilon_i^A=\epsilon_i^B$

对于k折交叉验证产生的k对测试错误率，先对每对的测试错误率求差，得到$\Delta_i=\epsilon_i^A-\epsilon_i^B$，若两个学习器性能相同，则差值均值应为0，因此对$\Delta_1,\Delta_2,…,\Delta_k$来对”学习器A和学习器B性能相同”的假设进行t验证,计算出$\Delta$的均值$\mu$和方差$\sigma^2$，在显著度$\alpha$下，若

小于临界值$t_{\alpha/2,k-1}$则假设不能被拒绝，认为两个学习器的性能没有显著差异；否则则认为两个学习器有显著差异，且平均错误率较小的学习器性能更好，其中$t_{\alpha/2,k-1}$为自由度为k-1的t分布上尾部累积分布为$\alpha/2$的临界值

进行有效的假设验证的一个重要的前提是测试错误率均为泛化错误率的独立 采样。但是通常情况下由于样本有限，在使用交叉验证等实验估计时，不同轮次的训练集都有不定程度的重叠，这使得测试错误率实际上并不独立，会导致过高估计假设成立的概率，为缓解这一问题，可采用”5$\times$2交叉验证”法。

就是做5次2折交叉验证，在每次2折交叉验证之前水机将数据打乱，使得5次交叉验证中的数据划分不重复。对于学习器A和B，第i次2折交叉验证将产生两对测试错误率，对他们分别求差，得到1折的差值$\Delta_i^1$和2折的差值$\Delta_i^2$，为了缓解测试错误率的非独立性，我们仅计算第1次2折交叉验证的两个结果的平均值，$\mu=0.5(\Delta_1^1+\Delta_1^2)$,但对每次2折实验的结果都计算其方差$\sigma_i^2=(\Delta_i^1-\frac{\Delta_i^1+\Delta_i^2}{2})^2+ (\Delta_i^2-\frac{\Delta_i^1+\Delta_i^2}{2})^2$,变量

服从自由度为5的t分布，其双边检验的临界值$t_{\alpha/2.5}$当，$\alpha=0.05$时，为2.5706，$\alpha=0.1$时，为2.0150

2.4.3 McNemar检验

对于二类问题，使用留出法不仅可估计出学习器A和B的测试错误率，还能获得两个学习器分类结果的差别

表: 两学习器分类差别列联表

假设两个学习器性能相同，则应有$e_{01}=e_{10}$，那么$\left|e_{01}-e_{10}\right|$应当服从正态分布，McNemar检验考察变量

因为$e_{01}+e_{10}非常小，需要考虑连续性校正，所以分子中有-1$

服从自由度为1的$\mathcal{X}^2$分布，即标准正态分布变量的平方，盖顶显著度$\alpha$，当以上变量值小于临界值$\mathcal{X}^2$时，不能拒绝假设，认为两个学习器性能没有显著区别，否则拒绝原假设，认为两个学习器有限制差别，并且平均错误率小的那个学习器性能更好。

2.4.4 Friedman检验与Nemenyi后续检验

交叉验证T检验和McNemar检验都是在一个数据集上比较两个算法的的性能，但是很多情况下会在一组数据集上对多个算法进行比较，当有多个算法参与比较时，种方法是在每个数据集上分别列出两两比较的结果，而在亮亮比较上可以用前面的方法；另一种方法更为直接，使用基于算法排序的Friedman检验

假设用$D_1,D_2,D_3,D_4$四个数据集对算法A、B、C进行比较，首先使用留出法或交叉验证法得到每个算法在每个数据集上的测试结果，然后在每个数据集上根据测试性能由好到坏排序，并赋值1,2,…;若算法的测试性能相同，则平分序值。列如得到如下的算法比较序值表

然后，使用Friedman检验来判断这些算法是否性能相同，若相同，则他们的平均序值应当相同，假设我们在N个数据集上比较k个算法，令$\mathcal{r}_i$表示第i个算法的平均序值，为了简化讨论，暂不考虑评分序值的情况，则$\mathcal{r}_i$的均值和方差分别为(k+1)/2和$(k^2-1)/12$.变量

在k和N都较大时，服从自由度为k-1的$\mathcal{X}^2$分布

然而上述的”原始Friedman检验”过于保守，因为原始检验需要k较大(例如k>30),现在通常使用变量

$\tau_F$服从自由度为k-1和(k-1)(N-1)的F分布

若”所有算法的性能相同”这个假设被拒绝，则说明算法的性能显著不同，这是需要进行”后续检验”(post-hoc test)来进一步区分各算法，常用的有Nemenyi后续检验

Nemenyi检验计算出平均序值差别的临界值域

其中qa为Tukey分布的临界值，可以通过算法个数k和$\alpha$这两个参数查表得到，若两个算法的平均序值之差超过了临界值域CD，则以相同的置信度拒绝”两个算法性能相同”这一假设

根据这一系列内容，可以得到上述例子的关系

求出$\tau_F$的值为24.429，查F检验的表，得到$\alpha=0.05$临界值为5.143，因为大于临界值，所以拒绝原假设，然后使用”后续检验”，根据$\alpha=0.05$和k=3得到CD值为1.657，根据平均序值的差，可以得到A与B、B与C这两对算法没有显著区别，但是A与C是有显著区别的

2.5偏差和方差

除了估计其泛化性能之外，还需要了解具有这样的泛化性能的本质，所以”偏差-方差分解”(bias-variance decomposition)是解释学习算法泛化性能的一种重要工具

偏差-方差分解试图对学习算法的期望泛化错误率进行拆解。对于测试样本x,令$y_D$为数据集中的标记，y为x的真实标记(因为有噪声，所以有可能$y_D\neq y$)$f(x;D)$为训练集D上学的模型f在x上的预测输出，以回归任务为例，

学习算法的期望预测为

使用样本数相同的不同训练集产生的方差为

噪声为

期望输出与实际标记的差别称为偏差(bias)，即

为了方便讨论，假定噪声期望为0，即$\epsilon^2=\mathbb{E}_D\left[(y_D-y)^2\right]=0$，通过简单的多项式展开合并，可以对算法的期望泛化误差进行分解

也就是说泛化误差可分解为偏差、方差和噪声之和

偏差度量了学习算法预测与真实结果的偏离程度，刻画了学习算法本身的拟合能力；方差度量了同样大小的训练集的变动所导致的学习性能的变化，即刻画了数据扰动所造成的影响；噪声则表达了在当前任务上任何学习算法能达到的期望泛化误差的下限，即刻画了学习问题本身的难度。

偏差-方差分解说明了泛化性能是由学习算法的能力、数据的充分性及学习任务本身的难度所共同决定的。为了获得好的泛化性能，则需要使偏差较小，即能够充分拟合数据，并且使得方差较小，即使得数据扰动影响小

一般而言，偏差和方差是有冲突的，称之为偏差-方差窘境(bias-variance dilemma)

一个训练模型主要有三个阶段，当训练不足时，因为学习器的拟合能力不强，所以训练数据的扰动不足以影响学习器，所以偏差主导泛化错误率，当训练的加深，学习器的拟合能力加强，训练数据的扰动渐渐被学习器所获取，方差主导泛化错误率，但随之逐步加深，学习器的拟合能力过强，导致一些数据局部性的特性被学习器学习到，导致发生过拟合，即把数据独有的特性当成了共同的特性。

习题

1. 数据集包含1000个样本，其中正例和反例都是500个，用留出法获取70%的训练数据和30%的测试数据，共有多少种划分方式？

留出法就是等比获取，也就是说在正例和反例都随机获取150个作为测试数据，不管训练数据，因为测试数据定好了，训练数据就是唯一得了。所以是$C_{150}^{500} C_{150}^{500}=\dbinom{500}{150}^2$个方式

1. 训练集包含100个样本，其中正反例各一半，假定学习算法所产生的模型是将新样本预测为训练样本较多的类别(训练样本数相同时进行随机猜测)，试给出10折交叉验证法和留一法分别对错误率进行估计所得的结果

k折交叉验证，因为分出的数据是尽量保持和原本数据集的分布一致性，所以都是正反例各一半，所以错误率是50%，

留一法，因为是选择一个数据样本当做测试集，所以根据训练样本反而因为测试集的那个样本的缺失，导致学习算法预测的是与测试集相反的结果，所以错误率是100%

1. 若学习器A的F1值比学习器B高，则A的BEP是否也比B高？

BEP是把训练得到的正例可能性将真实数据从高到低进行排序，在根据各个阀值进行阶段，得到的曲线中找到P和R一样的点。

F1则是根据混淆矩阵来的到P和R进行计算得到的。

根据https://blog.csdn.net/icefire_tyh/article/details/52065867 提供的例子可以看出BEP和F1并不想对等

可以看出因为根据学习器A和学习器B的排序，将真实数据从高到低进行排序，得到的顺序是一致的，所以BEP一致，但是按照F1计算，则可以得到

学习器A和B的混淆矩阵

学习器A和B的P、R、F1值

1. 试述正例率(TPR)，假正例率(FPR)，和查准率(P)、查全率(R)之间联系

TPR和R是一样的，都是预测的正确正例在所有真实正例的比例。

FPR是预测是正例但是实际是反例在所有真实反例的比例。P是预测正确的正例在所有预测为正例的比例，并没有实际关系

1. 试证明$AUC = 1 - \ell_{rank}$

首先解析AUC公式，AUC公式是根据求线下的面积得到的。

根据南瓜书(https://datawhalechina.github.io/pumpkin-book/#/chapter2/chapter2 )得知，书上的图为一个特例，即根据预测度进行排序后，根据每一个样本点是正例或者反例从而垂直或者水平前进；如果两个点都在阀值上，并且是一个正例和反例时，那则直线应该是斜方向前进的。不过无论是那种情况，都可以使用 $(上底+下底)高 \cfrac{1}{2}$来获取阀值间隔内的面积。

$\ell_{rank}=\frac{1}{m^+m^-}\sum_{x^+\in D^+}\sum_{x^-\in D^-}(I(f(x^+)<f(x^-))+ \frac{1}{2}I(f(x^+)=f(x^-)))$

根据南瓜书提供的图例，可以得到3个情况

绿线为垂直向上，表示1个正例被正确判断，在y轴上的投影为”$\frac{1}{m^+}$”

红线为水平向右，表示一个反例被错误判断成了正例，在x上的投影为”$\frac{1}{m^-}$”

蓝线为朝着右上方向的斜线，表示有a个正例和b个反例被判断成了正例，则在x轴的投影为”$b\frac{1}{m^-}$”,在y轴上的投影为”$a\frac{1}{m^+}$”

其中图像的每条折线的预测值相同，并且从0点开始逐渐较小。

将公式拆分开来，
$\sum_{x^+ \in D^+}$可以看做对所有正例的一个遍历，分解成

for $x_i^+ \;in \;D^+: \\
\qquad \frac{1}{m^+m^-}\sum_{x^-\in D^-}(I(f(x^+)<f(x^-))+ \frac{1}{2}I(f(x^+)=f(x^-)))\qquad 记为公式S$

其中$x_i^+$对应着每条绿线或者蓝线线段，遍历$x_i^+$可以看做遍历所有绿线和蓝线，S求出对应的绿线或蓝线和y轴组成的面积。

详细内容，先看绿线
对于每条绿线线段，因为是单一的绿线，所以没有预测点既有正例和反例的情况，所以去除后半段，只看反例的预测点大于正例的预测点的情况。
所以S简化为

其中每一个空格的长度为$\frac{1}{m^-}$高度为$\frac{1}{m^+}$,$\sum_{x^-\in D^-}I(f(x_i^+)<f(x^-))$为预测点大于正例的反例的个数，对应的就是空格的数量，乘以高和长，得到绿线左面的面积。

对于蓝线，公式S展开为

有两个部分，其中前面的一部分就是蓝线左边的空格的面积，而右边的为自身的格子的左边的面积的，即为自身格子的左边的三角形面积。

因为总面积为1，所以可以得到$AUC = 1 - \ell_{rank}$

1. 试述错误率和ROC曲线的关系

错误率为$\frac{FP+FN}{FP+FN+TP+TN}$

ROC有两个值为TPR和FPR，分别为$\frac{TP}{TP+FN}$，$\frac{FP}{FP+TN}$
每个点对应着一对TPR和FPR。

因为$m^+ = TP+FN,\;m^-=TN+FP$,所以用来代替，得到

所以每个点可以用错误率来代替。错误率最小的点为FPR为0，TPR为1的时候

代价敏感的错误率为$\frac{cost_{01}\times FN + cost_{10}\times FP}{m}$
所以代价敏感的错误率为

1. 证明任意一条ROC曲线都有一条对应代价曲线与之对应，反之亦然

以后补充

1. Min—max规范化和z-score规范化两种常见的规范化方式。这两者有何优缺点。
   令x为规范化前的取值，x’为规范化后的取值，其中$x_{min},x_{max}$为规范化前的最小值和最大值，$x_{min}’,x_{max}’$为规范化后的最小和最大取值，$\bar{x},\sigma_x$为规范化前的均值和标准差,则:

Min-max的规范化公式为

z-score的规范化公式为

Min-max规范化的优势是，每次新的数据如果在$x_{min}$ $x_{max}$的范围内，则不需要重新计算其他值，并且得到规范化的值都是大于0的;但是如果有极小(大)值，会有影响

z-score的优势是，极值对于规范化的影响不大，但是每次有新值时，都需要重新计算$\bar{x}$ $\sigma_x$

1. 试述$\mathcal{X}^2$检验过程

适用于定类和定类的检验。

以后补充

1. 试述Friedman检验中使用的两个方程区别

即$\tau_{\mathcal{X}^2}=\frac{12N}{k(k+1)}(\sum\limits_{i=1}^k\mathcal{r}_i^2-\frac{k(k+1)^2}{4})$

和$\tau_F=\frac{(N-1)\tau_{\mathcal{X}^2}}{N(k-1)-\tau_{\mathcal{X}^2}}$的区别

只知道$\tau_{\mathcal{X}^2}$需要满足k较大(k>30)的情况下，在k较小的情况下，会倾向于认为无显著区别。但是$\tau_F$并没有这个影响

详细的以后补充

Welcome to Hexo! This is your very first post. Check documentation for more info. If you get any problems when using Hexo, you can find the answer in troubleshooting or you can ask me on GitHub.

Quick Start

Create a new post

1

$ hexo new "My New Post"

More info: Writing

Run server

1

$ hexo server

More info: Server

Generate static files

1

$ hexo generate

More info: Generating

Deploy to remote sites

1

$ hexo deploy

More info: Deployment

支持向量机

分隔超平面为$w^Tx + b$
常数b类似于回归里的$w_0$

因为w为超平面的法向量，其中设距离为$\gamma$,$x_i 到x_0$的距离为$\gamma$,其中$x_0$为在平面的点，x为要求距离的点，

x的向量为$x_0 + \gamma \frac{w}{||w||}$ ,其中$\frac{w}{||w||}$为向量w的的单位向量，又有$w^Tx_0 + b = 0$与$w^Tw = ||w||^2$

$x-x_0 = \gamma \frac{\vec{w}}{||w||}$ 两边乘以$w^T$得

$w^T(x-x_0)= \gamma \frac{w^Tw}{||w||}$

$w^T(x-x_0)= \gamma \frac{||w||^2}{||w||}$

$w^Tx + b=\gamma||w||$

所以计算点A到平面的距离:值为$\frac{|w^TA + b| }{||w||}$

当$y_i$为1,-1时，$\gamma = \frac{2}{||w||}$ $s.t. w^Tx+b >= 1$

所以求$\gamma$的最大值，就是求$\frac{1}{2}||w||^2$ $s.t. w^Tx+b >= 1$的最小值

用拉格朗日乘子法来把约束条件加到求解方程里

$L(w,\alpha,b)=\frac{1}{2}||w||^2 + \sum\limits_{i=1}^m\alpha^i(1-y_i(w^Tx_i+b)$

令$\theta(w)=\mathop{max}\limits_{a_{i} >=0}L(w.\alpha,b)$

$\mathop{min}\limits_{w,b}\theta(w)=\mathop{min}\limits_{w,b}\mathop{max}\limits_{a_i >=0}L(w,\alpha,b)=p^*$

用对偶性，最大和最小互换，但是有kkt条件需要满足,才能保证$p^=d^$

$\mathop{max}\limits_{\alpha_i>=0}\mathop{mim}\limits_{w,b}L(w,\alpha,b)=d^*$

kkt条件为

min(f(x)),

$s.t. \\h_j(x)=0,j=1,2,..,p \\
g_k(x) <=0,k=1,2,…,q \\
x\in X\subset R^n
$

需要满足

1. $L(w,\alpha,b)$对x求导为0
2. $h_i(x)=0$
3. $\alpha * g_k(x) =0$

对w,或b求导得到
$\frac{\partial L}{\partial w}=0 得到 w= \sum\limits_{i=1}^n\alpha_iy_ix_i$

$\frac{\partial L}{\partial b} =0 得到 \sum\limits_{i=1}^n\alpha_iy_i=0$

把w带回L中得到，一个只有$\alpha$的方程,??为什么是减去

$\mathop{max}\limits_{\alpha}\sum\limits_{i=1}^n\alpha_i - \frac{1}{2}\sum\limits_{i,j=1}^{n}\alpha_i\alpha_jy_iy_jx_i^Tx_j$

$s.t. \alpha_i >=0,i=1,2,…n, \\
\sum\limits_{i=1}^n\alpha_iy_i=0$

SMO内容

$\alpha$的部分

将方程变形，乘以-1,变成一个最小的问题

$\mathop{min}\limits_{\alpha}\frac{1}{2}\sum\limits_{i=1}^{n}\sum\limits_{j=1}^n\alpha_i\alpha_jy_iy_jx_i^Tx_j -\sum\limits_{i=1}^n\alpha_i $

$s.t. \alpha_i >=0,i=1,2,…n, \\
\sum\limits_{i=1}^n\alpha_iy_i=0$

加入松弛变量C，保证可以有部分不是线性可分的，条件变为

$s.t. C>= \alpha_i >=0,i=1,2,…n, \\
\sum\limits_{i=1}^n\alpha_iy_i=0$

SMO为更新两个$\alpha$

则 $\alpha_1^{new}y_1 + \alpha_2^{new}y_2 = \alpha_1^{old}y_1 + \alpha_2^{old}y_2 = \zeta $

固定其他的值，只改变$\alpha_2^{new}的解$先确定它的上下限

$L<=\alpha_2^{new}<=H$

因为$C>=\alpha>=0$ 和$\alpha_1^{new}y_1 + \alpha_2^{new}y_2 = \alpha_1^{old}y_1 + \alpha_2^{old}y_2 = \zeta $

所以

当$y_1\ne y_2$时

$\alpha_1 - \alpha_2 = \zeta$

$\alpha_2 = \alpha_1 - \zeta$

得到

$C-\zeta>=\alpha_2 >= -\zeta$

结合$C>=\alpha>=0$,得到

$L=max(0,-\zeta),H=min(C,C-\zeta)$

同理，当$y_1=y_2$时，

$L=max(0,\zeta-C),H=min(C,\zeta)$

因此得到$\alpha_2^{new}$的上下界L和H为

$L=max(0,\alpha_2^{old}-\alpha_1^{old}),H=min(C,C+\alpha_2^{old} -\alpha_1^{old}),ify_1\ne y_2$

$L=max(0,\alpha_2^{old}+\alpha_1^{old}-C),H=min(C,\alpha_2^{old} +\alpha_1^{old}),ify_1= y_2$

固定除了$\alpha_1 \alpha_2的其他变量$

$w(\alpha_2) = \sum\limits_{i=1}^n\alpha_i - \frac{1}{2}\sum\limits_{i=1}^n\sum\limits_{i=1}^ny_iy_jx_i^Tx_j\alpha_i\alpha_j$

将$\alpha_1 \alpha_2$提出来得到

$w(\alpha_2) = \alpha_1 + \alpha_2 - \frac{1}{2}\alpha_1^2x_1^Tx_1 - \frac{1}{2}\alpha_2^2x_2^Tx_2 - y_1y_2\alpha_1\alpha_2x_1^Tx_2 - y_1\alpha_1v_1 - y_2\alpha_2v_2 + constant$

其中

定义$f(x_i) = \sum\limits_{j=1}^n\alpha_jy_jx_i^Tx_j + b$

$v_i =\sum\limits_{j=3}^n\alpha_jy_jx_i^Tx_j = f(x_i) -\sum\limits_{j=1}^2\alpha_jy_jx_i^Tx_j- b $

之后找到$\alpha_1和\alpha_2的关系$

因为$\sum\limits_{i=1}^n\alpha_iy_i=0$，所以将除了$\alpha_1y_1,\alpha_2y_2$的其他项看做常数-B

$\alpha_1y_1 + \alpha_2y_2 = B$，等式乘以$y_1$得到

$\alpha_1 = \gamma - s\alpha_2$ $，其中\gamma为By_1,s为y_1y_2$

带入公式，并且进行偏导

$\frac{\partial W(\alpha_2)}{\partial\alpha_2}=-s +1 +\gamma sx_1^Tx_1 - \alpha_2 x_1^Tx_1 - \alpha_2 x_2^Tx_2 - \gamma sx_1^Tx_2 + 2\alpha_2x_1^Tx_2 + y_2v_1 - y_2v_2 = 0$

导入$s=y_1y_2$得到

$\alpha_2^{new} = \frac{y_2(y_2-y_1+y_1\gamma(x_1^\intercal x_1 - x_1^\intercal x_2)+ v_1 -v_2)}{x_1^\intercal x_1 + x_2^\intercal x_2-2x_1^\intercal x_2}$

令$E_i = f(x_i)-y_i$

$\eta=x_1^Tx_1 + x_2^Tx_2 - 2x_1^Tx_2$

$E_i$为误差项，$\eta$为学习速率

已知$\gamma= \alpha_1^{old} + s\alpha_2^{old}$,和

$v_j=\sum\limits_{i=3}^n\alpha_iy_ix_j^Tx_i = f(x_j) -\sum\limits_{i=1}^2\alpha_iy_ix_j^Tx_i- b $$

简化$\alpha_2^{new}$

$\alpha_2^{new}=\alpha_2^{old} +\frac{y_2(E_1-E_2)}{\eta}$

加上约束，最终得到的解为

$\alpha_2^{new,clipped}=\begin{cases}
H,& \mbox{if }a_2^{new} >H \\
a_2^{new},& \mbox{if }L<=a_2^{new}<=H \\
L,& \mbox{if }a_2^{new} <L
\end{cases}$

又因为

$\alpha_1^{old}=\gamma -s\alpha_2^{old}$

$\alpha_1^{new} = \gamma -s\gamma_2^{new,clipped}$

消去$\gamma$得到

$\alpha_1^{new} = \alpha_1^{old} + y_1y_2(\alpha_2^{old} - \alpha_2^{new,clipped})$

b的部分

根据$y_1(w^Tx_1 + b)=1$两边乘以$y_1$，得到

$w^Tx_i+b =y_1$,又因为$w= \sum\limits_{i=1}^n\alpha_iy_ix_i$,得到

$\sum\limits_{i=1}^n\alpha_iy_ix_ix_1 +b=y1$单独提出$\alpha_1,alpha_2$，得

$b_1^{new} = y_1 -\sum\limits_{i=3}^n\alpha_iy_ix_i^Tx_1 - \alpha_1^{new}y_1x_1^Tx_1 - \alpha_2^{new}y_2x_2^Tx_1$

其中前两项为
$y_1 -\sum\limits_{i=3}^n\alpha_iy_ix_i^Tx_1=-E_1 + \alpha_1^{old}y_1x_1^Tx_1 + \alpha_2^{old}y_2x_2^Tx_1 + b^{old}$

整理得到

$b_1^{new}=b^{old} -E_1-y_1(\alpha_1^{new}- \alpha_1^{old})x_1^Tx_1 -y_2(\alpha_2^{new}-\alpha_2^{old})x_2^Tx_1$

同理得到$b_2^{new}$

$b_2^{new}=b^{old} -E_2-y_1(\alpha_1^{new}- \alpha_1^{old})x_1^Tx_2 -y_2(\alpha_2^{new}-\alpha_2^{old})x_2^Tx_2$

当$b_1,b_2都有效时，b^{new}=b_1^{new} = b_2^{new}$

所以，b的取值为

SMO计算梳理

1. 计算误差$E_i$

$E_i=f(x_i) - y_i = \sum\limits_{j=1}^n\alpha_jy_jx_i^Tx_j +b -y_i$

1. 计算上下限

1. 计算$\eta$

$\eta=x_1^Tx_1 + x_2^Tx_2 - 2x_1^Tx_2$

1. 更新$\alpha_j$

$\alpha_j^{new}=\alpha_j^{old} +\frac{y_j(E_1-E_j)}{\eta}$

1. 根据取值范围修剪$\alpha_j$

$\alpha_j^{new,clipped}=\begin{cases}
H,& \mbox{if }a_j^{new} >H \\
a_2^{new},& \mbox{if }L<=a_j^{new}<=H \\
L,& \mbox{if }a_j^{new} <L
\end{cases}$

1. 更新$\alpha_i$

$alpha_i^{new} = \alpha_i^{old} + y_jy_i(\alpha_j^{old} - \alpha_j^{new,clipped})$

1. 更新$b_i,b_j$

$b_i^{new}=b^{old} -E_i-y_i(\alpha_i^{new}- \alpha_i^{old})x_i^Tx_i -y_j(\alpha_j^{new}-\alpha_j^{old})x_j^Tx_i$

$b_j^{new}=b^{old} -E_j-y_i(\alpha_i^{new}- \alpha_i^{old})x_i^Tx_j -y_j(\alpha_j^{new}-\alpha_j^{old})x_j^Tx_j$

1. 根据bi,bj更新b

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28

import numpy as np
#SMO辅助函数
def loadDataset(file_name):
    """读取文件，得到两个list，
    data_mat 为2列m行，label_mat为m个元素"""
    data_mat=[];label_mat=[]
    fr = open(file_name)
    for line in fr.readlines():
        line_arr = line.strip().split('\t')
        data_mat.append([float(line_arr[0]),float(line_arr[1])])
        label_mat.append(float(line_arr[2]))
    return data_mat,label_mat

def selectJrand(i,m):
    """根据i，随机选取一个不是i的j"""
    import random
    j = i
    while(j==i):
        j = int(random.uniform(0,m))
    return j

def clipAlpha(aj,H,L):
    """修剪alpha"""
    if aj > H:
        aj = H
    elif aj < L:
        aj = L
    return aj

1
2

data_arr,label_arr = loadDataset('../../Downloads/machinelearninginaction/Ch06/testSet.txt')
len(data_arr)

100

1

set(label_arr)

{-1.0, 1.0}

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83

#流程为
#创建一个alpha向量，初始化为0
#当迭代次数小于最大迭代次数时(外循环)
#  对数据集的每一个数据向量(内循环)
#    如果该向量可以被优化:
#       随机选择另一个数据向量
#       同时优化这两个向量
#       如果两个向量都不能被优化，退出内循环
#  如果所有向量都没有被优化，增加迭代数目，继续下一次循环


def smo_simple(data_mat_in,class_labels,C,toler,max_iter):
    """数据集，数据标签，松弛常数，容错率，最大循环次数"""
    data_matrix = np.mat(data_mat_in);label_mat =np.mat(class_labels).transpose()
    b= 0;m,n = np.shape(data_matrix)
    alphas = np.mat(np.zeros((m,1)))
    iter = 0
    #外循环为总的循环数
    while(iter < max_iter):
        alpha_pairs_chanaged = 0
        #内循环为每个求E_i,判断是否要优化，能优化则随机找一个j，进行优化，优化成功则改变alpha_pairs_changed，
        #一次内循环后，如果有改变则改iter为0，继续进行，如果最大循环次数中都没有改变则跳出外循环，得到结果
        for i in range(m):
            #1.求E_i
            fx_i = float(np.multiply(alphas,label_mat).T * (data_matrix * data_matrix[i,:].T)) + b
            e_i = fx_i - float(label_mat[i])
            #根据kkt标准
            #alpha=0时，y_ig(x_i) >=1
            #0<alpha<C时，y_ig(x_i) = 1
            #alpha =C时，y_Ig(x_i) <= 1
            #如果不符合其中的一条，则进行优化
            if ((label_mat[i] * e_i < -toler) and (alphas[i] < C) or
               ((label_mat[i] * e_i > toler) and (alphas[i] > 0))):
                j = selectJrand(i,m)
                #求出了f(x_j)
                fx_j = float(np.multiply(alphas,label_mat).T * (data_matrix * data_matrix[j,:].T)) + b
                #fx_j = float(np.multiply(alphas,label_mat).T) * (data_matrix[j,:] * data_matrix.T).T + b #为原本方程的顺序
                e_j = fx_j - float(label_mat[j]) #得到E_j
                
                #2.求上下限
                alpha_i_old = alphas[i].copy()
                alpha_j_old = alphas[j].copy()
                if (label_mat[i] != label_mat[j]):
                    L = max(0,alphas[j] - alphas[i])
                    H = min(C,C+ alphas[j] - alphas[i])
                else:
                    L = max(0,alphas[j] + alphas[i] - C)
                    H = min(C, alphas[j] + alphas[i])
                if L==H: #上下限相等，C为0，或者都为C，或者一个为C，一个为0；C为0时，则alpha都是0，不成立；都为C时，没法增大也没法变小，一个为C，一个为0时，即alpha_i为C，并不会改变，还是C
                    #print('L==H')
                    continue
                #eta与原本的符号相反
                eta = 2.0 * data_matrix[i,:] * data_matrix[j,:].T - data_matrix[i,:] * data_matrix[i,:].T - data_matrix[j,:] * data_matrix[j,:].T
                if eta >= 0: #变化率为负数，并不符合，eta应为正数
                    #print('eta >=0')
                    continue
                #原本的加上，现在为减去
                alphas[j] -= label_mat[j] *(e_i - e_j)/eta
                alphas[j] = clipAlpha(alphas[j],H,L)
                if (abs(alphas[j] - alpha_j_old) <0.00001):
                    #print('j not moving enough') #变化太小，跳出本次循环，换其他j
                    continue
                alphas[i] += label_mat[i]*label_mat[j] *(alpha_j_old - alphas[j])
                
                b_1 = b - e_i - label_mat[i] * (alphas[i] - alpha_i_old) * data_matrix[i,:] * data_matrix[i,:].T -  label_mat[j] * (alphas[j] - alpha_j_old)* data_matrix[i,:] * data_matrix[j,:].T
                
                b_2 = b - e_j - label_mat[i] * (alphas[i] -alpha_i_old) * data_matrix[i,:] * data_matrix[j,:].T - label_mat[j] * (alphas[j] - alpha_j_old) * data_matrix[j,:] * data_matrix[j,:].T
                if (0 <alphas[i]) and (C > alphas[i]): 
                    b = b_1
                elif (0 < alphas[j]) and (C > alphas[j]):
                    b = b_2
                else:
                    b = (b_1 + b_2) /2.0
                alpha_pairs_chanaged += 1
                #print('iter: %d i:%d,pairs changed %d'%(iter,i,alpha_pairs_chanaged))
                #print('alpha_i_old: %f, alpha_i_new: %f,alpha_j_old: %f, alpha_j_new: %f' %(alpha_i_old,alphas[i],alpha_j_old,alphas[j]))
        print("fullSet, iter: %d , pairs changed %d" % (iter,alpha_pairs_chanaged))
        if (alpha_pairs_chanaged == 0):
            iter += 1
        else:
            iter =0
        print('iteration number: %d' % iter)
    return b,alphas

1

b,alphas = smo_simple(data_arr,label_arr,0.6,0.001,10)

fullSet, iter: 0 , pairs changed 7
iteration number: 0
fullSet, iter: 0 , pairs changed 4
iteration number: 0
fullSet, iter: 0 , pairs changed 0
iteration number: 1
fullSet, iter: 1 , pairs changed 6
iteration number: 0
fullSet, iter: 0 , pairs changed 1
iteration number: 0
fullSet, iter: 0 , pairs changed 2
iteration number: 0
fullSet, iter: 0 , pairs changed 1
iteration number: 0
fullSet, iter: 0 , pairs changed 1
iteration number: 0
fullSet, iter: 0 , pairs changed 0
iteration number: 1
fullSet, iter: 1 , pairs changed 0
iteration number: 2
fullSet, iter: 2 , pairs changed 0
iteration number: 3
fullSet, iter: 3 , pairs changed 1
iteration number: 0
fullSet, iter: 0 , pairs changed 1
iteration number: 0
fullSet, iter: 0 , pairs changed 0
iteration number: 1
fullSet, iter: 1 , pairs changed 0
iteration number: 2
fullSet, iter: 2 , pairs changed 0
iteration number: 3
fullSet, iter: 3 , pairs changed 1
iteration number: 0
fullSet, iter: 0 , pairs changed 0
iteration number: 1
fullSet, iter: 1 , pairs changed 1
iteration number: 0
fullSet, iter: 0 , pairs changed 0
iteration number: 1
fullSet, iter: 1 , pairs changed 2
iteration number: 0
fullSet, iter: 0 , pairs changed 0
iteration number: 1
fullSet, iter: 1 , pairs changed 2
iteration number: 0
fullSet, iter: 0 , pairs changed 0
iteration number: 1
fullSet, iter: 1 , pairs changed 0
iteration number: 2
fullSet, iter: 2 , pairs changed 0
iteration number: 3
fullSet, iter: 3 , pairs changed 1
iteration number: 0
fullSet, iter: 0 , pairs changed 1
iteration number: 0
fullSet, iter: 0 , pairs changed 1
iteration number: 0
fullSet, iter: 0 , pairs changed 1
iteration number: 0
fullSet, iter: 0 , pairs changed 0
iteration number: 1
fullSet, iter: 1 , pairs changed 0
iteration number: 2
fullSet, iter: 2 , pairs changed 0
iteration number: 3
fullSet, iter: 3 , pairs changed 0
iteration number: 4
fullSet, iter: 4 , pairs changed 1
iteration number: 0
fullSet, iter: 0 , pairs changed 0
iteration number: 1
fullSet, iter: 1 , pairs changed 1
iteration number: 0
fullSet, iter: 0 , pairs changed 0
iteration number: 1
fullSet, iter: 1 , pairs changed 0
iteration number: 2
fullSet, iter: 2 , pairs changed 1
iteration number: 0
fullSet, iter: 0 , pairs changed 0
iteration number: 1
fullSet, iter: 1 , pairs changed 0
iteration number: 2
fullSet, iter: 2 , pairs changed 0
iteration number: 3
fullSet, iter: 3 , pairs changed 2
iteration number: 0
fullSet, iter: 0 , pairs changed 0
iteration number: 1
fullSet, iter: 1 , pairs changed 0
iteration number: 2
fullSet, iter: 2 , pairs changed 0
iteration number: 3
fullSet, iter: 3 , pairs changed 0
iteration number: 4
fullSet, iter: 4 , pairs changed 0
iteration number: 5
fullSet, iter: 5 , pairs changed 0
iteration number: 6
fullSet, iter: 6 , pairs changed 0
iteration number: 7
fullSet, iter: 7 , pairs changed 1
iteration number: 0
fullSet, iter: 0 , pairs changed 0
iteration number: 1
fullSet, iter: 1 , pairs changed 0
iteration number: 2
fullSet, iter: 2 , pairs changed 0
iteration number: 3
fullSet, iter: 3 , pairs changed 0
iteration number: 4
fullSet, iter: 4 , pairs changed 0
iteration number: 5
fullSet, iter: 5 , pairs changed 0
iteration number: 6
fullSet, iter: 6 , pairs changed 0
iteration number: 7
fullSet, iter: 7 , pairs changed 0
iteration number: 8
fullSet, iter: 8 , pairs changed 0
iteration number: 9
fullSet, iter: 9 , pairs changed 0
iteration number: 10

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142

#边界上的样本对应的α_i=0或者α_i=C，在优化过程中很难变化，然而非边界样本0<α_i<C会随着对其他变量的优化会有大的变化
from numpy import *
class optStruct:
    def __init__(self,dataMatIn, classLabels, C, toler):  # Initialize the structure with the parameters 
        self.X = dataMatIn
        self.labelMat = classLabels
        self.C = C
        self.tol = toler
        self.m = shape(dataMatIn)[0]
        self.alphas = mat(zeros((self.m,1)))
        self.b = 0
        self.eCache = mat(zeros((self.m,2))) #第一项储存是否有效
        
def calcEk(oS, k):
    fXk = float(multiply(oS.alphas,oS.labelMat).T*(oS.X*oS.X[k,:].T)) + oS.b
    Ek = fXk - float(oS.labelMat[k])
    return Ek
        
def selectJ(i, oS, Ei):         #this is the second choice -heurstic, and calcs Ej
    maxK = -1; maxDeltaE = 0; Ej = 0
    oS.eCache[i] = [1,Ei]  #set valid #choose the alpha that gives the maximum delta E
    validEcacheList = nonzero(oS.eCache[:,0].A)[0]
    if (len(validEcacheList)) > 1:
        for k in validEcacheList:   #loop through valid Ecache values and find the one that maximizes delta E
            if k == i: continue #don't calc for i, waste of time
            Ek = calcEk(oS, k)
            deltaE = abs(Ei - Ek)
            if (deltaE > maxDeltaE):
                maxK = k; maxDeltaE = deltaE; Ej = Ek
        return maxK, Ej
    else:   #in this case (first time around) we don't have any valid eCache values
        j = selectJrand(i, oS.m)
        Ej = calcEk(oS, j)
    return j, Ej

def updateEk(oS, k):#after any alpha has changed update the new value in the cache
    Ek = calcEk(oS, k)
    oS.eCache[k] = [1,Ek]
        
def innerL(i, oS):
    Ei = calcEk(oS, i)
    #判断是否可以优化，不可以直接返回0，
    #可以优化则搜索j，第一次是随机的，并且添加0到临时e表里，因为只有一个数值，所以随机选择j
    #之后更新alpah-i j对，更新i,j的e值，返回结果
    #第二次进入内训，i=1,添加1到临时表，j因为有两个所以选j=0，
    #一遍循环后临时e表都有了，之后就能按照alpah不是0的来选择最佳的j
    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,Ej = selectJ(i, oS, Ei) #this has been changed from selectJrand
        alphaIold = oS.alphas[i].copy(); alphaJold = oS.alphas[j].copy();
        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 = 2.0 * oS.X[i,:]*oS.X[j,:].T - oS.X[i,:]*oS.X[i,:].T - oS.X[j,:]*oS.X[j,:].T
        if eta >= 0: 
            #print("eta>=0")
            return 0
        oS.alphas[j] -= oS.labelMat[j]*(Ei - Ej)/eta
        oS.alphas[j] = clipAlpha(oS.alphas[j],H,L)
        #因为alpha_j改变了，所以需要重新计算，为什么不是一起更新
        updateEk(oS, j) #added this for the Ecache
        if (abs(oS.alphas[j] - alphaJold) < 0.00001): 
            #print("j not moving enough")
            return 0
        oS.alphas[i] += oS.labelMat[j]*oS.labelMat[i]*(alphaJold - oS.alphas[j])#update i by the same amount as j
        updateEk(oS, i) #added this for the Ecache                    #the update is in the oppostie direction
        b1 = oS.b - Ei- oS.labelMat[i]*(oS.alphas[i]-alphaIold)*oS.X[i,:]*oS.X[i,:].T - oS.labelMat[j]*(oS.alphas[j]-alphaJold)*oS.X[i,:]*oS.X[j,:].T
        b2 = oS.b - Ej- oS.labelMat[i]*(oS.alphas[i]-alphaIold)*oS.X[i,:]*oS.X[j,:].T - oS.labelMat[j]*(oS.alphas[j]-alphaJold)*oS.X[j,:]*oS.X[j,:].T
        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 innerL1(i, oS): #一起更新临时e
    Ei = calcEk(oS, i)
    #判断是否可以优化，不可以直接返回0，
    #可以优化则搜索j，第一次是随机的，并且添加0到临时e表里，因为只有一个数值，所以随机选择j
    #之后更新alpah-i j对，更新i,j的e值，返回结果
    #第二次进入内训，i=1,添加1到临时表，j因为有两个所以选j=0，
    #一遍循环后临时e表都有了，之后就能按照alpah不是0的来选择最佳的j
    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,Ej = selectJ(i, oS, Ei) #this has been changed from selectJrand
        alphaIold = oS.alphas[i].copy(); alphaJold = oS.alphas[j].copy();
        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 = 2.0 * oS.X[i,:]*oS.X[j,:].T - oS.X[i,:]*oS.X[i,:].T - oS.X[j,:]*oS.X[j,:].T
        if eta >= 0: 
            #print("eta>=0"); 
            return 0
        oS.alphas[j] -= oS.labelMat[j]*(Ei - Ej)/eta
        oS.alphas[j] = clipAlpha(oS.alphas[j],H,L)
        if (abs(oS.alphas[j] - alphaJold) < 0.00001): 
            #print("j not moving enough")
            return 0
        oS.alphas[i] += oS.labelMat[j]*oS.labelMat[i]*(alphaJold - oS.alphas[j])#update i by the same amount as j
        updateEk(oS,j)
        updateEk(oS, i) #added this for the Ecache                    #the update is in the oppostie direction
        b1 = oS.b - Ei- oS.labelMat[i]*(oS.alphas[i]-alphaIold)*oS.X[i,:]*oS.X[i,:].T - oS.labelMat[j]*(oS.alphas[j]-alphaJold)*oS.X[i,:]*oS.X[j,:].T
        b2 = oS.b - Ej- oS.labelMat[i]*(oS.alphas[i]-alphaIold)*oS.X[i,:]*oS.X[j,:].T - oS.labelMat[j]*(oS.alphas[j]-alphaJold)*oS.X[j,:]*oS.X[j,:].T
        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):    #full Platt SMO
    oS = optStruct(mat(dataMatIn),mat(classLabels).transpose(),C,toler)
    iter = 0
    entireSet = True; alphaPairsChanged = 0
    #当循环数不满时并且是循环整个数据集或者不是循环整个数据集但是有alpha对改变时
    while (iter < maxIter) and ((alphaPairsChanged > 0) or (entireSet)):
        alphaPairsChanged = 0
        if entireSet:   #go over all
            for i in range(oS.m):        
                alphaPairsChanged += innerL(i,oS)
                #print("fullSet, iter: %d i:%d, pairs changed %d" % (iter,i,alphaPairsChanged))
            print("fullSet, iter: %d , pairs changed %d" % (iter,alphaPairsChanged))
            iter += 1
        else:#go over non-bound (railed) alphas
            nonBoundIs = 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))
            print("non-bound, iter: %d , pairs changed %d" % (iter,alphaPairsChanged))
            iter += 1
        #第一次循环为整个，循环后改变为循环non-bound，如果没有a对改变，则变成整个的数据集
        #在整个数据集循环后没有a对改变时，跳出循环，结束
        if entireSet: entireSet = False #toggle entire set loop
        elif (alphaPairsChanged == 0): entireSet = True  
        print("iteration number: %d" % iter)
    return oS.b,oS.alphas

1

data_arr, label_arr = loadDataset('../../Downloads/machinelearninginaction/Ch06/testSet.txt')

1

b, alphas = smoP(data_arr,label_arr,0.6,0.001,40)

fullSet, iter: 0 , pairs changed 6
iteration number: 1
non-bound, iter: 1 , pairs changed 0
iteration number: 2
fullSet, iter: 2 , pairs changed 0
iteration number: 3

$w = \sum\limits_{i=1}^m\alpha_iy_ix_i$

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

#计算w值
def calc_ws(alphas,data_arr,class_labels):
    x = np.mat(data_arr)
    label_mat = np.mat(class_labels).transpose()
    m,n = np.shape(x)
    w = np.zeros((n,1))
    for i in range(m):
        w += np.multiply(alphas[i] * label_mat[i], x[i,:].T)
    return w

def calc_ws(alphas,data_arr,class_labels):
    #向量来求解
    x = np.mat(data_arr);label_mat = np.mat(class_labels).T
    w = x.T * np.multiply(alphas,label_mat) 
    return w

1

calc_ws(alphas,data_arr,label_arr)

matrix([[ 0.65307162],
        [-0.17196128]])

1

calc_ws(alphas,data_arr,label_arr)

matrix([[ 0.65307162],
        [-0.17196128]])

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29

#绘制样本
def showClassifer(dataMat, w, b):
    #绘制样本点
    data_plus = []                                  #正样本
    data_minus = []                                 #负样本
    for i in range(len(dataMat)):
        if label_arr[i] > 0:
            data_plus.append(dataMat[i])
        else:
            data_minus.append(dataMat[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], s=30, alpha=0.7)   #正样本散点图
    plt.scatter(np.transpose(data_minus_np)[0], np.transpose(data_minus_np)[1], s=30, alpha=0.7) #负样本散点图
    #绘制直线
    x1 = max(dataMat)[0]
    x2 = min(dataMat)[0]
    a1, a2 = w
    b = float(b)
    a1 = float(a1[0])
    a2 = float(a2[0])
    y1, y2 = (-b- a1*x1)/a2, (-b - a1*x2)/a2
    plt.plot([x1, x2], [y1, y2])
    #找出支持向量点
    for i, alpha in enumerate(alphas):
        if abs(alpha) > 0:
            x, y = dataMat[i]
            plt.scatter([x], [y], s=150, c='none', alpha=0.7, linewidth=1.5, edgecolor='red')
    plt.show()

1
2

from matplotlib import pylab as plt
%matplotlib inline

c越高，说明越不能容忍出现误差,容易过拟合。C越小，容易欠拟合。C过大或过小，泛化能力变差

1
2
3

b, alphas = smo_simple(data_arr,label_arr,0.6,0.001,10)
ws = calc_ws(alphas,data_arr,label_arr)
showClassifer(data_arr,ws,b)

fullSet, iter: 0 , pairs changed 3
iteration number: 0
fullSet, iter: 0 , pairs changed 1
iteration number: 0
fullSet, iter: 0 , pairs changed 4
iteration number: 0
fullSet, iter: 0 , pairs changed 0
iteration number: 1
fullSet, iter: 1 , pairs changed 3
iteration number: 0
fullSet, iter: 0 , pairs changed 2
iteration number: 0
fullSet, iter: 0 , pairs changed 2
iteration number: 0
fullSet, iter: 0 , pairs changed 4
iteration number: 0
fullSet, iter: 0 , pairs changed 3
iteration number: 0
fullSet, iter: 0 , pairs changed 0
iteration number: 1
fullSet, iter: 1 , pairs changed 1
iteration number: 0
fullSet, iter: 0 , pairs changed 1
iteration number: 0
fullSet, iter: 0 , pairs changed 0
iteration number: 1
fullSet, iter: 1 , pairs changed 1
iteration number: 0
fullSet, iter: 0 , pairs changed 1
iteration number: 0
fullSet, iter: 0 , pairs changed 3
iteration number: 0
fullSet, iter: 0 , pairs changed 0
iteration number: 1
fullSet, iter: 1 , pairs changed 0
iteration number: 2
fullSet, iter: 2 , pairs changed 0
iteration number: 3
fullSet, iter: 3 , pairs changed 0
iteration number: 4
fullSet, iter: 4 , pairs changed 1
iteration number: 0
fullSet, iter: 0 , pairs changed 0
iteration number: 1
fullSet, iter: 1 , pairs changed 2
iteration number: 0
fullSet, iter: 0 , pairs changed 0
iteration number: 1
fullSet, iter: 1 , pairs changed 1
iteration number: 0
fullSet, iter: 0 , pairs changed 0
iteration number: 1
fullSet, iter: 1 , pairs changed 2
iteration number: 0
fullSet, iter: 0 , pairs changed 0
iteration number: 1
fullSet, iter: 1 , pairs changed 1
iteration number: 0
fullSet, iter: 0 , pairs changed 0
iteration number: 1
fullSet, iter: 1 , pairs changed 0
iteration number: 2
fullSet, iter: 2 , pairs changed 0
iteration number: 3
fullSet, iter: 3 , pairs changed 0
iteration number: 4
fullSet, iter: 4 , pairs changed 1
iteration number: 0
fullSet, iter: 0 , pairs changed 0
iteration number: 1
fullSet, iter: 1 , pairs changed 0
iteration number: 2
fullSet, iter: 2 , pairs changed 0
iteration number: 3
fullSet, iter: 3 , pairs changed 1
iteration number: 0
fullSet, iter: 0 , pairs changed 0
iteration number: 1
fullSet, iter: 1 , pairs changed 1
iteration number: 0
fullSet, iter: 0 , pairs changed 0
iteration number: 1
fullSet, iter: 1 , pairs changed 0
iteration number: 2
fullSet, iter: 2 , pairs changed 0
iteration number: 3
fullSet, iter: 3 , pairs changed 2
iteration number: 0
fullSet, iter: 0 , pairs changed 1
iteration number: 0
fullSet, iter: 0 , pairs changed 0
iteration number: 1
fullSet, iter: 1 , pairs changed 0
iteration number: 2
fullSet, iter: 2 , pairs changed 0
iteration number: 3
fullSet, iter: 3 , pairs changed 0
iteration number: 4
fullSet, iter: 4 , pairs changed 0
iteration number: 5
fullSet, iter: 5 , pairs changed 0
iteration number: 6
fullSet, iter: 6 , pairs changed 0
iteration number: 7
fullSet, iter: 7 , pairs changed 0
iteration number: 8
fullSet, iter: 8 , pairs changed 0
iteration number: 9
fullSet, iter: 9 , pairs changed 1
iteration number: 0
fullSet, iter: 0 , pairs changed 0
iteration number: 1
fullSet, iter: 1 , pairs changed 0
iteration number: 2
fullSet, iter: 2 , pairs changed 0
iteration number: 3
fullSet, iter: 3 , pairs changed 0
iteration number: 4
fullSet, iter: 4 , pairs changed 0
iteration number: 5
fullSet, iter: 5 , pairs changed 0
iteration number: 6
fullSet, iter: 6 , pairs changed 1
iteration number: 0
fullSet, iter: 0 , pairs changed 0
iteration number: 1
fullSet, iter: 1 , pairs changed 0
iteration number: 2
fullSet, iter: 2 , pairs changed 0
iteration number: 3
fullSet, iter: 3 , pairs changed 3
iteration number: 0
fullSet, iter: 0 , pairs changed 0
iteration number: 1
fullSet, iter: 1 , pairs changed 0
iteration number: 2
fullSet, iter: 2 , pairs changed 0
iteration number: 3
fullSet, iter: 3 , pairs changed 0
iteration number: 4
fullSet, iter: 4 , pairs changed 0
iteration number: 5
fullSet, iter: 5 , pairs changed 0
iteration number: 6
fullSet, iter: 6 , pairs changed 0
iteration number: 7
fullSet, iter: 7 , pairs changed 0
iteration number: 8
fullSet, iter: 8 , pairs changed 0
iteration number: 9
fullSet, iter: 9 , pairs changed 0
iteration number: 10

1
2
3

b, alphas = smoP(data_arr,label_arr,0.6,0.001,40)
ws = calc_ws(alphas,data_arr,label_arr)
showClassifer(data_arr,ws,b)

fullSet, iter: 0 , pairs changed 6
iteration number: 1
non-bound, iter: 1 , pairs changed 0
iteration number: 2
fullSet, iter: 2 , pairs changed 0
iteration number: 3

1
2
3

b, alphas = smoP(data_arr,label_arr,0.01,0.001,40)
ws = calc_ws(alphas,data_arr,label_arr)
showClassifer(data_arr,ws,b)

fullSet, iter: 0 , pairs changed 10
iteration number: 1
non-bound, iter: 1 , pairs changed 0
iteration number: 2
fullSet, iter: 2 , pairs changed 0
iteration number: 3

1
2
3

b, alphas = smoP(data_arr,label_arr,55,0.001,40)
ws = calc_ws(alphas,data_arr,label_arr)
showClassifer(data_arr,ws,b)

fullSet, iter: 0 , pairs changed 8
iteration number: 1
non-bound, iter: 1 , pairs changed 0
iteration number: 2
fullSet, iter: 2 , pairs changed 0
iteration number: 3

核函数
径向基核函数
$k(x,y) = exp(\frac{-||x -y||^2}{2\sigma^2})$

自变量为向量，可以计算自变量相对于(0,0)或者其他的向量距离

$\sigma$为到达率(reach)或者函数值跌落到0的速度参数

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34

#转换核函数
def kernelTrans(X,A,kTup):
    """X为整个数据集，A为其中一行，计算A行和整个数据集的核函数内积
        如果是线性，的则直接得到x * xi^T
    """
    m,n = shape(X)
    K = np.mat(np.zeros((m,1)))
    if kTup[0] == 'lin':
        K = X * A.T #数据集i列的内积 (m,n) * (n,1) => (m,1)
    elif kTup[0] == 'rbf':
        for j in range(m): #读取第j行的值,x 
            deltaRow = X[j,:]  -A #j行减去i行的差，x-y (1,n)
            K[j] = deltaRow * deltaRow.T #得到L2反数
        K = exp(K / (-1 * kTup[1] ** 2)) #各值计算得到和函数的值,k有m行
    else:
        raise NameError('Houston we have a problem --That Kernel is not recognized')
    return K


class optStructk:
    def __init__(self,dataMatIn, classLabels, C, toler, kTup):  # Initialize the structure with the parameters 
        self.X = dataMatIn
        self.labelMat = classLabels
        self.C = C
        self.tol = toler
        self.m = shape(dataMatIn)[0]
        self.alphas = mat(zeros((self.m,1)))
        self.b = 0
        self.eCache = mat(zeros((self.m,2))) #first column is valid flag
        #K为i,j的内积
        self.K = mat(zeros((self.m,self.m)))
        for i in range(self.m):
            #得到第i列的值即y
            self.K[:,i] = kernelTrans(self.X, self.X[i,:], kTup)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62

#修改用到核函数的辅助函数
def calcEk(oS, k):
    fXk = float(multiply(oS.alphas,oS.labelMat).T*oS.K[:,k] + oS.b)
    Ek = fXk - float(oS.labelMat[k])
    return Ek
    
def innerLk(i, oS):
    Ei = calcEk(oS, i)
    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,Ej = selectJ(i, oS, Ei) #this has been changed from selectJrand
        alphaIold = oS.alphas[i].copy(); alphaJold = oS.alphas[j].copy();
        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 = 2.0 * oS.K[i,j] - oS.K[i,i] - oS.K[j,j] #changed for kernel
        if eta >= 0: 
            #print("eta>=0"); 
            return 0
        oS.alphas[j] -= oS.labelMat[j]*(Ei - Ej)/eta
        oS.alphas[j] = clipAlpha(oS.alphas[j],H,L)
        updateEk(oS, j) #added this for the Ecache
        if (abs(oS.alphas[j] - alphaJold) < 0.00001):
            #print("j not moving enough"); 
            return 0
        oS.alphas[i] += oS.labelMat[j]*oS.labelMat[i]*(alphaJold - oS.alphas[j])#update i by the same amount as j
        updateEk(oS, i) #added this for the Ecache                    #the update is in the oppostie direction
        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 smoPk(dataMatIn, classLabels, C, toler, maxIter,kTup=('lin', 0)):    #full Platt SMO
    oS = optStructk(mat(dataMatIn),mat(classLabels).transpose(),C,toler, kTup)
    iter = 0
    entireSet = True; alphaPairsChanged = 0
    while (iter < maxIter) and ((alphaPairsChanged > 0) or (entireSet)):
        alphaPairsChanged = 0
        if entireSet:   #go over all
            for i in range(oS.m):        
                alphaPairsChanged += innerLk(i,oS)
                #print("fullSet, iter: %d i:%d, pairs changed %d" % (iter,i,alphaPairsChanged))
            print("fullSet, iter: %d , pairs changed %d" % (iter,alphaPairsChanged))
            iter += 1
        else:#go over non-bound (railed) alphas
            nonBoundIs = nonzero((oS.alphas.A > 0) * (oS.alphas.A < C))[0]
            for i in nonBoundIs:
                alphaPairsChanged += innerLk(i,oS)
                #print("non-bound, iter: %d i:%d, pairs changed %d" % (iter,i,alphaPairsChanged))
            print("non-bound, iter: %d , pairs changed %d" % (iter,alphaPairsChanged))
            iter += 1
        if entireSet: entireSet = False #toggle entire set loop
        elif (alphaPairsChanged == 0): entireSet = True  
        print("iteration number: %d" % iter)
    return oS.b,oS.alphas

$f(x) = \sum\limits_{i=1}^m\alpha_i^*y_iK(x\cdot x_i) + b^*$

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

def testRbf(k1=1.3):
    dataArr,labelArr = loadDataset('../../Downloads/machinelearninginaction/Ch06/testSetRBF.txt')
    b,alphas = smoPk(dataArr, labelArr, 200, 0.0001, 10000, ('rbf', k1)) #C=200 important
    datMat=mat(dataArr); labelMat = mat(labelArr).transpose()
    #得到支持向量的索引
    svInd=nonzero(alphas.A>0)[0] #[0]是因为alphas是一个mat，[1]是都是0的array
    sVs=datMat[svInd] #get matrix of only support vectors
    labelSV = labelMat[svInd];
    print("there are %d Support Vectors" % shape(sVs)[0])
    m,n = shape(datMat)
    errorCount = 0
    for i in range(m):
        #除了支持向量的a都是0，只需要计算支持向量的核函数
        kernelEval = kernelTrans(sVs,datMat[i,:],('rbf', k1)) #计算支持向量与数据的和函数，(shape(sVs)[0],1)     
        predict=kernelEval.T * multiply(labelSV,alphas[svInd]) + b #(1,i) * (i,1) =>(1,1)
        if sign(predict)!=sign(labelArr[i]): errorCount += 1
    print("the training error rate is: %f" % (float(errorCount)/m))
    dataArr,labelArr = loadDataset('../../Downloads/machinelearninginaction/Ch06/testSetRBF2.txt')
    errorCount = 0
    datMat=mat(dataArr); labelMat = mat(labelArr).transpose()
    m,n = shape(datMat)
    for i in range(m):
        kernelEval = kernelTrans(sVs,datMat[i,:],('rbf', k1))
        predict=kernelEval.T * multiply(labelSV,alphas[svInd]) + b
        if sign(predict)!=sign(labelArr[i]): errorCount += 1    
    print("the test error rate is: %f" % (float(errorCount)/m) )

1

testRbf()

fullSet, iter: 0 , pairs changed 30
iteration number: 1
non-bound, iter: 1 , pairs changed 5
iteration number: 2
non-bound, iter: 2 , pairs changed 2
iteration number: 3
non-bound, iter: 3 , pairs changed 0
iteration number: 4
fullSet, iter: 4 , pairs changed 0
iteration number: 5
there are 29 Support Vectors
the training error rate is: 0.130000
the test error rate is: 0.150000

$\sigma$ 如果太小，会得到很多支持向量，因为各支持向量的影响会变小，所以需要更多支持向量，但是容易过拟合

$\sigma $如果过大，则支持向量变小，容易欠拟合

手写识别问题回顾

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48

def img2vector(filename):
    returnVect = zeros((1,1024))
    fr = open(filename)
    for i in range(32):
        lineStr = fr.readline()
        for j in range(32):
            returnVect[0,32*i+j] = int(lineStr[j])
    return returnVect

def loadImages(dirName):
    from os import listdir
    hwLabels = []
    trainingFileList = listdir(dirName)           #load the training set
    m = len(trainingFileList)
    trainingMat = zeros((m,1024))
    for i in range(m):
        fileNameStr = trainingFileList[i]
        fileStr = fileNameStr.split('.')[0]     #take off .txt
        classNumStr = int(fileStr.split('_')[0])
        if classNumStr == 9: hwLabels.append(-1)
        else: hwLabels.append(1)
        trainingMat[i,:] = img2vector('%s/%s' % (dirName, fileNameStr))
    return trainingMat, hwLabels    

def testDigits(kTup=('rbf', 10)):
    dataArr,labelArr = loadImages('../../Downloads/machinelearninginaction/Ch06/digits/trainingDigits')
    b,alphas = smoPk(dataArr, labelArr, 200, 0.0001, 10000, kTup)
    datMat=mat(dataArr); labelMat = mat(labelArr).transpose()
    svInd=nonzero(alphas.A>0)[0]
    sVs=datMat[svInd] 
    labelSV = labelMat[svInd];
    print("there are %d Support Vectors" % shape(sVs)[0])
    m,n = shape(datMat)
    errorCount = 0
    for i in range(m):
        kernelEval = kernelTrans(sVs,datMat[i,:],kTup)
        predict=kernelEval.T * multiply(labelSV,alphas[svInd]) + b
        if sign(predict)!=sign(labelArr[i]): errorCount += 1
    print("the training error rate is: %f" % (float(errorCount)/m))
    dataArr,labelArr = loadImages('../../Downloads/machinelearninginaction/Ch06/digits/testDigits')
    errorCount = 0
    datMat=mat(dataArr); labelMat = mat(labelArr).transpose()
    m,n = shape(datMat)
    for i in range(m):
        kernelEval = kernelTrans(sVs,datMat[i,:],kTup)
        predict=kernelEval.T * multiply(labelSV,alphas[svInd]) + b
        if sign(predict)!=sign(labelArr[i]): errorCount += 1    
    print("the test error rate is: %f" % (float(errorCount)/m))

1

testDigits()

fullSet, iter: 0 , pairs changed 104
iteration number: 1
non-bound, iter: 1 , pairs changed 11
iteration number: 2
non-bound, iter: 2 , pairs changed 0
iteration number: 3
fullSet, iter: 3 , pairs changed 0
iteration number: 4
there are 115 Support Vectors
the training error rate is: 0.000000
the test error rate is: 0.016129

axis

轴用来为超过一维的数组定义的属性，二维数据拥有两个轴：第0轴沿着行的垂直往下，第1
轴沿着列的方向水平延伸。 逐行即把每一行延伸下来的一列当做一组
逐列即把每一列延伸下来的一行当做一组

1
2
3

df = pd.DataFrame(np.arange(12).reshape((3,4)),index =['a','b','c'],
                  columns =['one','two','three','four'] )
print(df)

1
2
3
4
5

#+RESULTS:
:    one  two  three  four
: a    0    1      2     3
: b    4    5      6     7
: c    8    9     10    11

sum等操作

1
2
3

print(df.sum()) #默认是0，逐行向下，所以是沿着列计算的
print()
print(df.sum(1))#axis = 1的操作，沿着列向右做计算的，所以是各行的计算

drop等的计算

1
2

print(df.drop('four',axis = 1)) #沿着列丢弃four
print(df.drop('c',axis = 0)) #沿着行丢弃c

1
2
3
4
5
6
7
8

#+RESULTS:
:    one  two  three
: a    0    1      2
: b    4    5      6
: c    8    9     10
:    one  two  three  four
: a    0    1      2     3
: b    4    5      6     7

plus拼接

1
2
3
4
5
6
7
8
9
10

df1 = df.reset_index().copy()
df2 =df.copy()

#merge
merged = df1.merge(df.reset_index(), on='index')
print(merged)  #横向连接,即axis = 1默认是内连接，
merged_index = df1.merge(df,right_index = True,left_on='index')
print(merged_index) #即df1的index列和df的index连接
merged_all_Index = df2.merge(df,right_index = True,left_index = True )
print(merged_all_Index)

1
2
3
4
5
6

#numpy的concatenate axis = 1即横向连接，axis=0即竖向连接,默认axis=0
arr = np.arange(12).reshape(3,4)
print(arr)
print(np.concatenate([arr,arr]))
print('axis = 1')
print(np.concatenate([arr,arr],axis = 1))

1
2
3
4

#pd.concat axis=0默认，即竖向连接，axis=1是横向，默认是outer连接
print(pd.concat([df,df]))
print('axis =1 ')
print(pd.concat([df,df],axis =1))

add/div等操作

1
2
3
4
5
6

#axis =0 的情况，即逐行的以组series
#因为df是 3*4的，所以要一个长度3的数组
print(df.add(np.arange(3),axis=0))

#axis = 1的情况， 即逐列的一组数组， 默认axis = 1
print(df.add(np.arange(4)))

1
2
3
4
5
6
7
8
9

#+RESULTS:
:    one  two  three  four
: a    0    1      2     3
: b    5    6      7     8
: c   10   11     12    13
:    one  two  three  four
: a    0    2      4     6
: b    4    6      8    10
: c    8   10     12    14

聚合groupby

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

#groupby 默认按照axis=0 也可以axis=1聚合
df = pd.DataFrame({'key1' : ['a', 'a', 'b', 'b', 'a'],
                   'key2' : ['one', 'two', 'one', 'two', 'one'],
                   'data1' : np.random.randn(5),
                   'data2' : np.random.randn(5)})
print(df)
#axis = 0即对于列key，划分数据，列1里，0，1，4为a,其他列同理；按照行分来
for a, b in df.groupby('key1',axis =0):
    print(a)
    print(b)
    print()

print('-----axis=1的情况----------')
#axis = 1 即对着列分割，每一行按照key进行分割，整体看来是按照列 分割
#object就是第一和第二列，float64即后面两列
for a, b in df.groupby(df.dtypes, axis =1):
    print(a)
    print(b)
    print()

chapter 1

基本import

1
2

import numpy as np
import pandas as pd

Chapter 2

Chapter 3: python

Chapter 4: numpy

创建ndarray

arange

1
2
3
4
5
6
7
8
9
10
11
12
13

import pandas as pd
pd.options.display.notebook_repr_html = False
import numpy as np
data1 =[6,7,8,0,1]
arr1 = np.array(data1)
arr1
print data1 * 2 #list 会变成重复数组
print arr1 * 2 # array会传播进行运算
#print data1 + 2 #会报错
print arr1 + 2
print arr1.ndim #维度
print arr1.shape
print arr1.reshape(5,1).shape

1
2
3
4
5
6
7
8

#+RESULTS:
    : [6, 7, 8, 0, 1, 6, 7, 8, 0, 1]
    : [12 14 16  0  2]
    : [ 8  9 10  2  3]
    : 1
    : (5,)
    : (5, 1)

创建函数

a = np.identity(3)
print(a)
a = np.ones_like(a)
print(a)

1
2
3
4
5
6
7

    #+RESULTS:
: [[ 1.  0.  0.]
:  [ 0.  1.  0.]
:  [ 0.  0.  1.]]
: [[ 1.  1.  1.]
:  [ 1.  1.  1.]
:  [ 1.  1.  1.]]

数据类型转换astype

1. 参数为数据类型

查看数据类型

1. array的自带属性dtype

索引和切片

切片是原始数据的视图，而不是复制，修改会反映到原数组上 注意点</span> [切片是原始数据的视图而不是复制修改会反映到原数组上]

新的表示方法

1
2
3
4
5
6
7
8

a = np.arange(15).reshape(3,5)
print(a)
print
print(a[2][3]) #原本的表示方法
print
print(a[2,3])  #numpy新增的表示方法
print
print a[[2,1]] #花式索引

1
2
3
4
5
6
7
8
9
10

[[ 0  1  2  3  4]
 [ 5  6  7  8  9]
 [10 11 12 13 14]]

13

13

[[10 11 12 13 14]
 [ 5  6  7  8  9]]

布尔型索引

1. 否定符号 != / -

   1 2 3 4 5 6 7 8 9 10 11

   names = np.array(["Bob","Joe", 'Will', 'Bob', "Will",'Joe','Joe']) data = np.random.randn(7,4) print(names) print print(data) print print(data[names == "Bob"]) print print(data[names != "Bob"]) print print(data[-(names == 'Bob')]) #numpy 新的表示方法

   1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25

   ['Bob' 'Joe' 'Will' 'Bob' 'Will' 'Joe' 'Joe'] [[ -1.06602469e+00 1.70960938e+00 -5.37441602e-01 -1.26864579e+00] [ 6.83392505e-01 -3.65908339e-01 -1.07189314e+00 1.78582204e+00] [ -1.06516810e+00 1.74430728e+00 5.88376963e-01 1.59791560e-01] [ -1.96586484e-01 1.14215100e-04 -7.00807814e-01 -1.75005941e-01] [ 9.63746755e-01 -4.51637875e-01 -1.78610575e+00 6.75211511e-01] [ -6.90968770e-01 5.80293483e-01 7.34267063e-01 4.58773330e-01] [ -1.36539167e-01 1.90051999e+00 -1.13351041e+00 7.74997645e-01]] [[ -1.06602469e+00 1.70960938e+00 -5.37441602e-01 -1.26864579e+00] [ -1.96586484e-01 1.14215100e-04 -7.00807814e-01 -1.75005941e-01]] [[ 0.6833925 -0.36590834 -1.07189314 1.78582204] [-1.0651681 1.74430728 0.58837696 0.15979156] [ 0.96374676 -0.45163788 -1.78610575 0.67521151] [-0.69096877 0.58029348 0.73426706 0.45877333] [-0.13653917 1.90051999 -1.13351041 0.77499765]] [[ 0.6833925 -0.36590834 -1.07189314 1.78582204] [-1.0651681 1.74430728 0.58837696 0.15979156] [ 0.96374676 -0.45163788 -1.78610575 0.67521151] [-0.69096877 0.58029348 0.73426706 0.45877333] [-0.13653917 1.90051999 -1.13351041 0.77499765]] /Users/Yueyec/.pyenv/versions/anaconda2-4.3.0/lib/python2.7/site-packages/ipykernel/__main__.py:11: DeprecationWarning: numpy boolean negative, the `-` operator, is deprecated, use the `~` operator or the logical_not function instead.

花式索引

1. 是复制，不是视图

2. 索引参数是一个list 里面是row index，可重复

3. 参数是二阶list，得到的会是对应的单个元素 注意点</span>

   1 2 3 4 5 6 7 8 9

   arr = np.arange(32).reshape((8,4)) print(arr) print print(arr[[1,5,7,2],[0,3,1,2]]) #会得到arr[1,0],arr[5,3],arr[7,1],arr[2,2]的一阶array print #第一种方法 print(arr[[1,5,7,2]][:,[0,3,1,2]]) print print(arr[np.ix_([1,5,7,2],[0,3,1,2])]) #np.ix_将两个一维整数数组转换成一个用于区域的索引器

   1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

   [[ 0 1 2 3] [ 4 5 6 7] [ 8 9 10 11] [12 13 14 15] [16 17 18 19] [20 21 22 23] [24 25 26 27] [28 29 30 31]] [ 4 23 29 10] [[ 4 7 5 6] [20 23 21 22] [28 31 29 30] [ 8 11 9 10]] [[ 4 7 5 6] [20 23 21 22] [28 31 29 30] [ 8 11 9 10]]

数组转置和轴对换

1. 转置为视图，非复制,但是不改变实例本身（self）的结构

2. T属性

   1. 轴颠倒

3. transpose()方法 没明白</span>

   1. 参数为由轴编号组成的元组

      1 2 3 4

      arr = np.arange(16).reshape((2,2,4)) print(arr) print "---------------------" print(arr.transpose((1,0,2)))

      1 2 3 4 5 6 7 8 9 10 11

      [[[ 0 1 2 3] [ 4 5 6 7]] [[ 8 9 10 11] [12 13 14 15]]] --------------------- [[[ 0 1 2 3] [ 8 9 10 11]] [[ 4 5 6 7] [12 13 14 15]]]

4. swapaxes函数 没明白</span>

通用函数

常用一元函数列表

常用二元函数列表

利用数组进行数据处理

假设在一组值（网格型）上计算函数sqrt(x^2^ + y^2^)

1
2
3
4
5
6

a = np.arange(5)
b = np.arange(5,9)
xs,ys = np.meshgrid(a, b)
print(xs)  #xs为a在轴0上广播len(b)次，即4
print '----------'
print(ys)  #ys为b.T在轴1上广播Len(a)次，即5

1
2
3
4
5
6
7
8
9

[[0 1 2 3 4]
 [0 1 2 3 4]
 [0 1 2 3 4]
 [0 1 2 3 4]]
----------
[[5 5 5 5 5]
 [6 6 6 6 6]
 [7 7 7 7 7]
 [8 8 8 8 8]]

1
2
3
4
5
6
7
8

points = np.arange(-5,5,0.01) #1000个间隔相等的点
xs,ys = np.meshgrid(points, points)
import matplotlib.pyplot as plt
z = np.sqrt(xs ** 2+ ys ** 2)
plt.imshow(z,cmap = plt.cm.gray)
plt.colorbar()
plt.title("Image plot of $\sqrt{x^2 + y^2}$ for a grid of values")
plt.show()

将条件逻辑表述为数组运算

numpy.where函数是三元表达式x if condition else y的矢量化版本

1. 假设这三个变量

   1 2 3 4

   import numpy as np xarr = np.array([1.1, 1.2, 1.3, 1.4, 1.5]) yarr = np.array([2.1, 2.2, 2.3, 2.4, 2.5]) cond = np.array([True, False, True, True, False])

2. 假设根据cond的值来选取xarr与yarr的值，当cond为True，选取xarr,否者yarr

   1. 常规做法，对于大数组处理速度不快，无法作用于多维数组

      1 2	result = [x if c else y for x, y, c in zip(xarr, yarr, cond)] print(result)

      1	[1.1000000000000001, 2.2000000000000002, 1.3, 1.3999999999999999, 2.5]

   2. 用where处理

      1 2	result = np.where(cond, xarr, yarr) print(result)

      1	[ 1.1 2.2 1.3 1.4 2.5]

3. 复杂的例子

   1. 常规的做法

      1 2 3 4 5 6 7 8 9 10 11 12 13

      cond1 = cond cond2 = np.array([True,False, False, True, True]) result = [] for i in range(len(cond1)): if cond1[i] and cond2[i]: result.append(0) elif cond1[i]: result.append(1) elif cond2[i]: result.append(2) else: result.append(3) print(result)

      1	[0, 3, 1, 0, 2]

   2. 用where

      1 2 3 4

      result = np.where(cond1 &cond2, 0, np.where(cond1, 1, np.where(cond2, 2, 3))) print(result)

      1	[0 3 1 0 2]

   3. 用“布尔值在计算过程中可以被当做0或1处理”，用算术运算（感觉像黑魔法）
      :有问题？:

      1 2	result = 1* (cond1 - cond2) + 2* (cond2 & - cond1) +3 * - (cond1 | cond2) print(result)

      1	[0 3 1 0 3]

数学和统计方法

sum, mean, std既可以当做实例方法，也可以当做顶级numpy函数使用

1
2
3
4
5
6
7

arr = np.random.randn(5,4) #正态分布的数据

print arr.mean()
print "-----------------------"
print np.mean(arr)
print "-----------------------"
print arr.sum()

1
2
3
4
5

0.0748211510483
-----------------------
0.0748211510483
-----------------------
1.49642302097

mean 与 sum这类函数可以接受一个axis参数（用于计算该轴上的统计值），最终结构是少一维度的数组 axis和想象中相反</span> [mean-与-sum这类函数可以接受一个axis参数用于计算该轴上的统计值最终结构是少一维度的数组]

1
2
3
4
5
6
7

print arr
print "---------------------"
print arr.mean(axis = 1) #在轴1中传播,即在列中进行计算
print "---------------------"
print arr.sum(0) #轴0
print "---------------------"
print arr.sum(1) #轴1

1
2
3
4
5
6
7
8
9
10
11

[[ 0.79860446 -0.67843272  2.58233563  0.33419308]
 [ 0.22848276 -1.42518949 -1.04403965 -0.3835565 ]
 [-0.34612276  1.39126673  0.29382674  0.53407341]
 [ 1.21116404 -1.05950897  0.18198447 -0.99575602]
 [ 0.81490393 -2.09189052  0.27434089  0.87574352]]
---------------------
[ 0.75917511 -0.65607572  0.46826103 -0.16552912 -0.03172555]
---------------------
[ 2.70703243 -3.86375497  2.28844807  0.36469749]
---------------------
[ 3.03670044 -2.62430288  1.87304412 -0.66211648 -0.12690218]

cumsum 和 cumprod之类不聚合，产生一个由中间结果组成的数组

1
2
3
4
5
6
7

import numpy as np
arr = np.array([[0,1,2],[3,4,5],[6,7,8]])
print(arr.cumsum(0))
print "-------------------"
print(arr.cumsum(1))
print "-------------------"
print(arr.cumprod(1))  #积累相乘

1
2
3
4
5
6
7
8
9
10
11

[[ 0  1  2]
 [ 3  5  7]
 [ 9 12 15]]
-------------------
[[ 0  1  3]
 [ 3  7 12]
 [ 6 13 21]]
-------------------
[[  0   0   0]
 [  3  12  60]
 [  6  42 336]]

基本数组统计方法

用于布尔型数组的方法

1. 可以用sum来求True的计数

   1 2	arr = np.random.randn(100) print((arr> 0).sum())

   1

   52

2. any all 来测试数组是否存在一个及以上的True,是否所有值为True

   1 2 3

   bools = np.array([False,False,True,False]) print(bools.any()) print(bools.all())

   1 2	True False

排序

1. numpy数组自带sort方法

   1 2 3 4 5

   arr = np.random.randn(8) print (arr) print "---------------------------------------------------" arr.sort() print (arr)

   1 2 3 4 5

   [ 1.8328002 0.41074931 -0.60400869 1.53465964 -0.60704556 -0.16515945 -1.02657211 -0.91581907] --------------------------------------------------- [-1.02657211 -0.91581907 -0.60704556 -0.60400869 -0.16515945 0.41074931 1.53465964 1.8328002 ]

2. 多维数组中任意轴上进行排序，只需参数线路轴编号

   1 2 3 4 5

   arr = np.random.randn(5,3) print arr print "-----------------------------" arr.sort(1) print arr

   1 2 3 4 5 6 7 8 9 10 11

   [[-0.51009631 -1.66867772 -0.22805943] [-0.66434219 0.38498122 -1.23044479] [-0.39114451 0.01949044 -1.25717288] [-0.00524041 0.42130187 -2.21689239] [ 0.95358967 -0.94485412 -0.84467746]] ----------------------------- [[-1.66867772 -0.51009631 -0.22805943] [-1.23044479 -0.66434219 0.38498122] [-1.25717288 -0.39114451 0.01949044] [-2.21689239 -0.00524041 0.42130187] [-0.94485412 -0.84467746 0.95358967]]

3. 顶级方法np.sort不会改变原数组，返回已经排序的副本

唯一化以及其他的集合逻辑

1. np.unique 返回数组已排序的唯一值

   1 2 3 4 5 6

   names = np.array(["Bob", 'Joe','Will','Bob','Will','Joe','Joe']) print np.unique(names) ints = np.array([3,3,3,2,2,1,1,4,4]) print np.unique(ints) print "-------------" print sorted(set(names)) #纯python代码，与之等价

   1 2 3 4

   ['Bob' 'Joe' 'Will'] [1 2 3 4] ------------- ['Bob', 'Joe', 'Will']

2. np.in1d 返回一个数组的值是否在另一个数组中，返回一个布尔型数组

   1 2 3 4

   values = np.array([6,0,0,3,2,5,6]) a = np.in1d(values, [2,3,6]) print a print values[a]

   1 2	[ True False False True True False True] [6 3 2 6]

3. 数组集合运算

   | 函数 | 描述 |
   |—————————|———————————————————————————————|
   | unique(x) | 计算x中的唯一元素，并返回有序的结果 |
   | intersect1d(x,y) | 计算x和y的公共元素，并返回有序结果 |
   | union1d(x,y) | 计算x和y的并集，并返回结果 |
   | in1d(x,y) | 得到一个“x的元素是否包含于y”的布尔型数组 |
   | setdiff1d(x,y) | 集合的差，机元素在x中且不在y中 |
   | setxor1d(x,y) | 集合的对称差，即存在于一个数组中但不同是存在两个数组中的元素 |
   | | |

用于数组的文件输入输出

1. 以数组为二进制格式保存到磁盘

   1. np.save, np.load

      1 2 3

      arr = np.arange(10) np.save('some_array',arr) #不加后缀名，自动添加.npy print np.load("some_array.npy")

      1	[0 1 2 3 4 5 6 7 8 9]

   2. np.savez
      可以将doge数组保存到一个压缩文件中，将数组以关键字参数的形式传入

      1 2 3 4 5

      arr1 = np.arange(10) arr2 = np.arange(1,20,2) np.savez("array_archive.npz",a = arr1, b = arr2) arch = np.load("array_archive.npz") print arch["b"]

2. 存取文本文件

   1. np.loadtxt(“xx”, deliniter = “,”) #第二个参数为分隔符

   2. np.genfromtxt() 面向的是结构化数组和缺失数据处理

线性代数

1. dot 方法

2. linalg

   1. 有一组标准的矩阵分解运算以及求逆和行列式 之类的东西

      1 2 3 4 5 6 7 8 9 10 11 12 13 14

      import numpy as np from numpy.linalg import inv, qr from numpy.random import randn x = randn(5,5) mat = x.T.dot(x) print "inv(mat)----------------" print inv(mat) print "mat.dot(inv(mat))---------------------------" print mat.dot(inv(mat)) q, r = qr(mat) print "q---------------------------" print q print "r---------------------------" print r

      1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34

      inv(mat)---------------- [[ 5.24583468 14.98884503 9.71335736 -7.77651755 -22.8396013 ] [ 14.98884503 46.52511417 28.81764808 -22.01122276 -69.70405132] [ 9.71335736 28.81764808 18.67257302 -14.00870932 -44.03024424] [ -7.77651755 -22.01122276 -14.00870932 12.70341866 32.9416411 ] [ -22.8396013 -69.70405132 -44.03024424 32.9416411 106.1381638 ]] mat.dot(inv(mat))--------------------------- [[ 1.00000000e+00 -1.27763158e-14 -3.16524540e-15 5.52810747e-15 -1.55455587e-14] [ -8.09523246e-15 1.00000000e+00 -5.76325857e-15 -1.66274227e-15 2.52694539e-14] [ 1.00097108e-15 2.71639602e-16 1.00000000e+00 -9.81680658e-15 3.54119545e-15] [ -6.64139148e-16 1.57971687e-14 3.29027055e-15 1.00000000e+00 -2.99822967e-14] [ -3.04200410e-15 -8.84058697e-15 1.92099465e-15 -8.57631485e-15 1.00000000e+00]] q--------------------------- [[-0.91326408 0.01491589 -0.30364189 0.21681123 -0.16285184] [-0.08041303 -0.74311076 0.43827531 -0.04710955 -0.49700662] [ 0.33046062 -0.31708316 -0.83157252 0.01337763 -0.31394621] [-0.21986762 -0.15288857 -0.13263187 -0.92494467 0.23488181] [-0.04400528 -0.56890005 -0.08131595 0.30833745 0.75679059]] r--------------------------- [[ -1.40092662e+01 -1.26845899e+00 5.84338504e+00 -3.35403461e+00 -3.83026371e-01] [ 0.00000000e+00 -2.13935952e+00 -2.42468320e+00 -6.51800989e-01 -2.21389665e+00] [ 0.00000000e+00 0.00000000e+00 -2.44765043e+00 -3.82142992e-01 -8.97542800e-01] [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 -4.11638083e-01 1.30663382e-01] [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 7.13024012e-03]]

   2. linalg的一些常用函数

      | 函数 | 描述 |
      |———-|————————————————————————————————————————————————-|
      | diag | 以一维数组的形式返回方阵的对角线（或非对角线）元素，或将一维数组转换成为方阵（非对角线元素为0） |
      | dot | 矩阵乘法 |
      | trace | 计算对角线元素的和 |
      | det | 计算矩阵行列式 |
      | eig | 计算方阵的本征值和本征向量 |
      | inv | 计算方阵的逆 |
      | pinv | 计算矩阵的Moore-Penrose违逆 |
      | qr | 计算QR分解 |
      | svd | 计算奇异值分解（SVD) |
      | solve | 解线性方程Ax= b,其中A为一个方阵 |
      | lstsq | 计算Ax = b的最小二乘解 |

随机数生成

1. numpy.random模块对内置random进行了补充

   1 2 3

   samples = np.random.normal(size = (4,4)) print samples

   1 2 3 4

   [[-0.2577954 -1.84871926 -1.73680389 0.75279352] [-0.39841909 -1.57836749 0.58627939 -0.02390301] [ 0.26546489 1.05297095 0.19336693 0.90791299] [-0.74260739 -1.06363081 0.27326372 -0.00851723]]

2. 部分numpy.random函数

   | 函数 | 描述 |
   |——————-|————————————————————————————————|
   | seed | 确定随机数生成器的种子 |
   | permutation | 返回一个序列的随机排列或返回一个随机排列的范围 |
   | shuffle | 对一个序列就地随机排列 |
   | rand | 产生均匀分布的样本值 |
   | randint | 从给定的上下限范围内随机选取整数,size参数（int) |
   | randn | 产生正态分布（平均值为0，标准差为1）的样本值，类似于MATLAB接口 |
   | binomial | 产生二项分布的样本值 |
   | normal | 产生正态（高斯）分布的样本值 |
   | beta | 产生Beta分布的样本值 |
   | chisquare | 产生卡方分布的样本值 |
   | gamma | 产生Gamma分布的样本值 |
   | uniform | 产生在[0,1)中均匀分布的样本值 |
   | | |

3. 范例：随机漫步

   1. 单个随机漫步

      1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22

      import numpy as np import random #传统python代码 position = 0 walk = [position] steps = 1000 for i in xrange(steps): step = 1 if random.randint(0,1) else -1 #新的写法 position += step walk.append(position) #numpy.random方法 nsteps = 1000 draws = np.random.randint(0,2,size = nsteps) steps = np.where(draws >0, 1,-1) walk = steps.cumsum() print "walk.min()----------" print walk.min() print "walk.max()------------" print walk.max()

      1 2 3 4

      walk.min()---------- -14 walk.max()------------ 26

   2. 一次性模拟多个随机漫步

      1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

      nwalks = 5000 nsteps = 1000 draws = np.random.randint(0,2, size = (nwalks, nsteps)) steps = np.where(draws >0, 1, -1) walks = steps.cumsum(1) print "walks.max()----------" print walks.max() print "walks.min()----------" print walks.min() #选出超过30或-30的walk hits30 = (np.abs(walks) >= 30).any(1) #any是在轴1上广播，即计算各个行中有没有绝对值30以上的 print "到达30或-30的数量" print hits30.sum() #到达30或-30的数量 #利用这个数据来选出那些行，并调用argmax在轴1上获取穿越时间 crossing_times = (np.abs(walks[hits30]) >= 30).argmax(1) print "crossing_times.mean()-------" print crossing_times.mean() #其他分布方式得到漫步数据 steps = np.random.normal(loc = 0, scale = 0.25, size = (nwalks,nsteps))

      1 2 3 4 5 6 7 8

      walks.max()---------- 119 walks.min()---------- -107 到达30或-30的数量 3329 crossing_times.mean()------- 497.677080204

Chapter 5: pandas入门

数据结构

Series

1. 类似于一维数组的对象，由一组数组和一组与之对应的数据标签（即索引）组成

   1 2 3

   import pandas as pd obj = pd.Series([4,7,-5,3]) print obj

   1 2 3 4 5

   0 4 1 7 2 -5 3 3 dtype: int64

2. values 和 index属性

   1 2	print obj.values print obj.index

   1 2	[ 4 7 -5 3] RangeIndex(start=0, stop=4, step=1)

3. 自定义index

   1 2 3

   obj2 = pd.Series([4,7,-5,3], index = ['d', 'b', 'a', 'c']) print obj2 print obj2.index

   1 2 3 4 5 6

   d 4 b 7 a -5 c 3 dtype: int64 Index([u'd', u'b', u'a', u'c'], dtype='object')

4. 通过索引来选取Series中的值

   1 2 3 4

   print obj2['a'] obj2['d'] = 6 print obj2[['c','a','d']]

   1 2 3 4 5

   -5 c 3 a -5 d 6 dtype: int64

   **

   1 2 3

   print obj2[obj2 > 0] print obj2 * 2 print np.exp(obj2)

   1 2 3 4 5 6 7 8 9 10 11 12 13 14

   d 6 b 7 c 3 dtype: int64 d 12 b 14 a -10 c 6 dtype: int64 d 403.428793 b 1096.633158 a 0.006738 c 20.085537 dtype: float64

5. 通过字典直接建立Series

   1 2 3

   sdata = {"Ohio" : 35000, 'Texas': 71000, "Oregon": 16000, "Utah": 5000} obj3 = pd.Series(sdata) print obj3

   1 2 3 4 5

   Ohio 35000 Oregon 16000 Texas 71000 Utah 5000 dtype: int64

6. 传入字典，并有index参数，则用index来排序,从字典中寻找index对应的值，没有的话，NaN来代替

   1 2 3

   states = ['California', "Ohio", 'Oregon','Texas'] obj4 = pd.Series(sdata,index = states) print obj4

   1 2 3 4 5

   California NaN Ohio 35000.0 Oregon 16000.0 Texas 71000.0 dtype: float64

7. pd.isnull, pd.notnull 寻找缺失数据

   1 2 3

   print pd.isnull(obj4) print print pd.notnull(obj4)

   1 2 3 4 5 6 7 8 9 10 11

   California True Ohio False Oregon False Texas False dtype: bool California False Ohio True Oregon True Texas True dtype: bool

8. Series 算术运算中会自动对齐不同的数据 注意点</span>

   1 2 3 4 5

   print obj3 print print obj4 print print obj3 + obj4

   1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18

   Ohio 35000 Oregon 16000 Texas 71000 Utah 5000 dtype: int64 California NaN Ohio 35000.0 Oregon 16000.0 Texas 71000.0 dtype: float64 California NaN Ohio 70000.0 Oregon 32000.0 Texas 142000.0 Utah NaN dtype: float64

9. Series 本身 与 index(索引） 有name属性

   1 2 3

   obj4.name = 'population' obj4.index.name = "state" print obj4

   1 2 3 4 5 6

   state California NaN Ohio 35000.0 Oregon 16000.0 Texas 71000.0 Name: population, dtype: float64

10. 修改 Series的索引

    1 2 3 4

    print obj print obj.index = ['Bob', "Steve", 'Jeff', "Ryan"] print obj

    1 2 3 4 5 6 7 8 9 10 11

    0 4 1 7 2 -5 3 3 dtype: int64 Bob 4 Steve 7 Jeff -5 Ryan 3 dtype: int64