本策略是为了验证SVM对于大盘涨跌的预测是否有效，相比于纯随机策略，是否有明显的提高。
SVM模型用06~14年的数据训练，16年1月~12月的数据用来回测，这样是为了避免因为在模型中投入了现阶段的数据导致的过拟合。
从结果来看，在16年的小范围下跌中，该模型表现还凑合吧……
import numpy as np
import pandas as pd
from CAL.PyCAL import Date
from CAL.PyCAL import Calendar
from CAL.PyCAL import BizDayConvention
from sklearn import svm
start = ''
# 回测起始时间
# 回测结束时间
benchmark = 'HS300'
# 策略参考标准
universe = set_universe('HS300')
# 证券池，支持股票和基金
capital_base = 100000
# 起始资金
freq = 'd'
# 策略类型，'d'表示日间策略使用日dw线回测，'m'表示日内策略使用分钟线回测
fields = ['tradeDate','closeIndex', 'highestIndex','lowestIndex', 'turnoverVol','CHG','CHGPct']
stock = '000300'
#tradeDate是交易日、closeIndex是收盘指数、highestIndex是当日最大指数，lowestIndex是当日最小指数，CHG是涨跌
index_raw = DataAPI.MktIdxdGet(ticker=stock,beginDate=u&&,endDate=u&&,field=fields,pandas=&1&)
#获取日到日，上一行代码设定的所有索引的相关信息。
index_date = index_raw.set_index('tradeDate')
index_date = index_date.dropna()
index_date['max_difference'] = index_date['highestIndex'] - index_date['lowestIndex']
index_date['max_of_30day'] = None
index_date['min_of_30day'] = None
index_date['max_difference_of_30day'] = None
index_date['closeIndex_after30days'] = None
#预设需要处理的值为None，方便之后直接用dropna函数去掉无效数据
for i in xrange(len(index_date)-30):
#对数据进行处理
index_date['max_of_30day'][i+30] = max(index_date['highestIndex'][i:i+30])
#找出前30天最大值。
index_date['min_of_30day'][i+30] = min(index_date['lowestIndex'][i:i+30])
#找出前30天最小值
index_date['max_difference_of_30day'][i+30] = max(index_date['max_difference'][i:i+30])
#找出前30天最大日波动
index_date['closeIndex_after30days'][i]=index_date['closeIndex'][i+30]
#找出30天后的收盘价。
index_date = index_date.dropna()
#去掉前30个和后30个无效的数据。
lables_raw = index_date['closeIndex_after30days'] #提取出需要预测的数据
lables = index_date['closeIndex_after30days'] & index_date['closeIndex'] #为分类处理数据，判断30天后的收盘价是否大于今日收盘价
lables_ud = lables.replace({True:'up',False:'down'}) #方便他人阅读，将True和False改为up和down，意味着30天后收盘价涨了还是跌了
features = index_date.drop(['closeIndex_after30days'],axis = 1) #在特征值中去掉我们要预测的数据。
在未调参之前，我们先获取一次准确率：
from sklearn import cross_validation
from sklearn import preprocessing
scaler = preprocessing.StandardScaler().fit(features)
features_scaler = scaler.transform(features)
#上面两行代码用来标准化数据
X_train,X_test, y_train, y_test = cross_validation.train_test_split(features_scaler, lables, test_size = 0.2, random_state = 0)
clf_svm = svm.SVC()
#使用SVM分类器来判断涨跌
clf_svm.fit(X_train, y_train)
print &预测准确率为：%0.2f& % (clf_svm.score(X_test, y_test))
然后调C值，这里我是先让C在1~100的range跑，然后100~200……到300~400的时候发现不会进一步提高了。其实可以直接从1~1000跑一次，很快就能绘画出整个变动的图，然而我电脑渣带不动。
i_list = []
score_list = []
for i in range(300,400,1):
clf_svm = svm.SVC(C = i)
#使用SVM分类器来判断涨跌
clf_svm.fit(X_train, y_train)
i_list.append(i)
score_list.append(clf_svm.score(X_test, y_test))
score_list_df =
pd.DataFrame({'i_list':i_list,'score_list':score_list})
score_list_df.plot(x='i_list' ,y='score_list',title='score change with c')
然后是gamma值，和C值调参上也是同理。
i_list = []
score_list = []
for i in range(100,200,1):
clf_svm = svm.SVC(C=350 ,
gamma = i)
#使用SVM分类器来判断涨跌
clf_svm.fit(X_train, y_train)
i_list.append(i)
score_list.append(clf_svm.score(X_test, y_test))
score_list_df =
pd.DataFrame({'gamma_list':i_list,'score_list':score_list})
score_list_df.plot(x='gamma_list' ,y='score_list',title='score change with gamma')
虽说没什么卵用……还是假吧意思的比对一下不同核函数下的准确率吧。理所当然的是默认的高斯核表现最好。
i_list = []
score_list = []
['linear', 'rbf','sigmoid']
for i in kernels :
clf_svm = svm.SVC(C=350 , gamma = 1.8 , kernel = i )
clf_svm.fit(X_train, y_train)
i_list.append(i)
score_list.append(clf_svm.score(X_test, y_test))
score_list_df =
pd.DataFrame({'kernels':i_list,'score_list':score_list})
score_list_df.plot(x='kernels' ,y='score_list',title='score change with kernels',kind='bar')
from sklearn.grid_search import GridSearchCV
from sklearn.cross_validation import ShuffleSplit
i = range(100,200,1)
params = {'C':range(300,400,1),'gamma':[x /100. for x in range(100,200,1)]}
# X_train,X_test, y_train, y_test = cross_validation.train_test_split(features_scaler, lables, test_size = 0.2, random_state = 0)
clf_svm = svm.SVC()
# cv_sets = ShuffleSplit(X_train.shape[0], n_iter = 10, test_size = 0.20, random_state = 0)
grid = GridSearchCV(clf_svm, params )
grid = grid.fit(X_train, y_train)
print grid.best_estimator_
然后在最优解的基础上再次计算一次准确率
from sklearn import cross_validation
from sklearn import preprocessing
scaler = preprocessing.StandardScaler().fit(features)
features_scaler = scaler.transform(features)
#上面两行代码用来标准化数据
X_train,X_test, y_train, y_test = cross_validation.train_test_split(features_scaler, lables, test_size = 0.2, random_state = 0)
clf_svm = svm.SVC(C = 300,gamma = 1.03)
#使用SVM分类器来判断涨跌
clf_svm.fit(X_train, y_train)
print &预测准确率为：%0.2f& % (clf_svm.score(X_test, y_test))
为了判断模型是否稳健，我们让训练集合处于变化中，然后观察随着训练集合的变化，准确率的波动范围图。这里采取的是数据每10个变化一次。
发现最低没有低于过0.72的准确率，波动较大在0.14左右，模型稳健度一般。
num_list = []
score_list = []
for i in xrange((len(features_scaler)-1000)/10):
num_now = len(features_scaler)%10 + 10*i +1000
X_train,X_test, y_train, y_test = cross_validation.train_test_split(features_scaler[:num_now], lables[:num_now], test_size = 0.2, random_state = 0)
clf_svm = svm.SVC(C=350,gamma = 1.8)
#使用SVM分类器来判断涨跌
clf_svm.fit(X_train, y_train)
num_list.append(num_now)
score_list.append(clf_svm.score(X_test, y_test))
score_list_df =
pd.DataFrame({'sets_num':num_list,'accuracy':score_list})
score_list_df.plot(x='sets_num' ,y='accuracy',title='Accuracy with sets')
接下来是比对用的空白组，纯随机策略（不控制风险，只是随机买，1.20倍卖出）
import random
start = ''
# 回测起始时间
# 回测结束时间
benchmark = 'HS300'
# 策略参考标准
universe = set_universe('HS300')
# 证券池，支持股票和基金
capital_base = 100000
# 起始资金
freq = 'd'
# 策略类型，'d'表示日间策略使用日线回测，'m'表示日内策略使用分钟线回测
refresh_rate = 1
# 调仓频率，表示执行handle_data的时间间隔，若freq = 'd'时间间隔的单位为交易日，若freq = 'm'时间间隔为分钟
def initialize(account):
# 初始化虚拟账户状态
features_list = []
def handle_data(account):
random.shuffle(account.universe)
# 随机化股票池一遍随机策略
for stock in account.universe:
# 股票是股票池中的股票，并且优矿帮你自动剔除了当天停牌退市的股票
p = account.reference_price[stock]
# 股票前一天的收盘价
cost = account.security_cost.get(stock)
# 股票的平均持仓成本
if not cost:
# 判断当前没有买入该股票
order_pct_to(stock, 0.10)
# 将满足条件的股票买入，总价值占虚拟帐户的10%
elif cost and p &= cost * 1.20:
# 卖出条件，当p这个价格涨幅到买入价的1.20倍；
order_to(stock, 0)
# 将满足条件的股票卖到剩余0股，即全部卖出
然后是纯随机策略基础上，只增加一个预测盘指的涨跌，如果预测涨，则随机买入，否则不买。和纯随机策略比，的确好了一丢丢。
import random
start = ''
# 回测起始时间
# 回测结束时间
benchmark = 'HS300'
# 策略参考标准
universe = set_universe('HS300')
# 证券池，支持股票和基金
capital_base = 100000
# 起始资金
freq = 'd'
# 策略类型，'d'表示日间策略使用日线回测，'m'表示日内策略使用分钟线回测
refresh_rate = 1
# 调仓频率，表示执行handle_data的时间间隔，若freq = 'd'时间间隔的单位为交易日，若freq = 'm'时间间隔为分钟
stock = '000300' #预测的指数，沪深300指数。和策略参考池一致。
fields = ['tradeDate','closeIndex', 'highestIndex','lowestIndex', 'turnoverVol','CHG','CHGPct']
#tradeDate是交易日、closeIndex是收盘指数、highestIndex是当日最大指数，lowestIndex是当日最小指数，CHG是涨跌
def initialize(account):
# 初始化虚拟账户状态
features_list = []
def handle_data(account):
# 生成买入列表
last_date = account.previous_date.strftime(&%Y-%m-%d&) #获取上一个交易日日期并格式化
begin_date = pd.date_range(end=last_date,periods=60)[0] #获取60日之前的交易日日期
begin_date = begin_date.strftime(&%Y-%m-%d&) #格式化这个日期
to_class = DataAPI.MktIdxdGet(ticker='000300',beginDate=begin_date,endDate=last_date,field=fields,pandas=&1&)
to_class = to_class.dropna()
to_class = to_class[-30:] #获取我们要的30天的指数信息
to_class_date = to_class.set_index('tradeDate')
to_class_date['max_difference'] = to_class_date['highestIndex'] - to_class_date['lowestIndex']
to_class_date_max_of_30day = max(to_class_date['highestIndex'])
#找出前30天最大值。
to_class_date_min_of_30day = min(to_class_date['lowestIndex'])
#找出前30天最小值
to_class_date_max_difference_of_30day = max(to_class_date['max_difference'])
#找出前30天最大日波动
features_for_predict = to_class_date[-1:]
features_for_predict['max_of_30day'] = to_class_date_max_of_30day
features_for_predict['min_of_30day'] = to_class_date_min_of_30day
features_for_predict['max_difference_of_30day'] = to_class_date_max_difference_of_30day
features_fp_scaler = scaler.transform(features_for_predict)
predict_up = clf_svm.predict(features_fp_scaler)
#预测30天后的收盘是涨还是跌。
random.shuffle(account.universe)
for stock in account.universe:
# 股票是股票池中的股票，并且优矿帮你自动剔除了当天停牌退市的股票
p = account.reference_price[stock]
# 股票前一天的收盘价
cost = account.security_cost.get(stock)
# 股票的平均持仓成本
if predict_up and not cost:
# 判断当前没有买入该股票
order_pct_to(stock, 0.10)
# 将满足条件的股票买入，总价值占虚拟帐户的10%
elif cost and p &= cost * 1.20:
# 卖出条件，当p这个价格涨幅到买入价的1.20倍；
order_to(stock, 0)
# 将满足条件的股票卖到剩余0股，即全部卖出
本文已收录于以下专栏：
相关文章推荐
Python股票历史数据预处理（一）
在进行量化投资交易编程时，我们需要股票历史数据作为分析依据，下面介绍如何通过Python获取股票历史数据并且将结果存为DataFrame格式。处理后的股票历史数...
网上有很多关于LibSVM的教程，LibSVM提供了多平台下的可执行文件，在win平台下可以直接通过命令行执行，可以利用提供的可执行文件方便的进行离线的训练和预测，这里不再赘述。这篇文章要做的是把SV...
from __future__ import division
from __future__ import print_function
import numpy as...
GoogLeNet, 2014 年 ILSVRC 挑战赛冠军,这个 model 证明了一件事：用更多的卷积，更深的层次可以得到更好的结构。结构拓宽卷积网络的宽度和深度，其中将稀疏矩阵合并成稠密矩阵的方...
机器学习算法与Python实践之（二）支持向量机（SVM）初级http://blog.csdn.net/zouxy09        机器学习算法与Python实践这个系列...
SVM 支持向量机
原理就不赘述了，其余的文章有讲过。SVM是一种十分优秀的分类算法，使用SVM也能给股票进行一定程度上的预测。
因为是分类算法，因此不像ARIMA一样预测的是时序...
SVM 支持向量机
如何使用libsvm进行回归预测
/thread-.html
(出处: MATLAB技术论坛)
强烈建议您看这个...
机器学习算法与Python实践这个系列主要是参考《机器学习实战》这本书。因为自己想学习Python，然后也想对一些机器学习算法加深下了解，所以就想通过Python来实现几个比较常用的机器学习算法。恰好...
他的最新文章
讲师：王哲涵
讲师：韦玮
您举报文章：
举报原因：
原文地址：
原因补充：
(最多只允许输入30个字)Scikit-Learn-1.2：支持向量机
来源：博客园
支持向量机是一种用于分类，回归和离群检测的一种监督学习的方法。
支持向量机的优势有：
高纬度有效
维度数大于样点数
1.2.1 分类
SVM有三种分类模型：SVC(C-Support Vector Classification.),NuSVC(Nu-Support Vector Classification.),LinearSVC(Linear Support Vector Classification).三者有一定差异。
SVC和NuSVC这两种方法相类似。但是他们接受有些不同的参数集。有不同的数学公式。LinearSVC是另一种支持向量分类的实施方式，它支持线性核函数的输入。他也相比SVC和NuSVC而言少了一些参量，如support_
这三种分类模型，都是输入两个阵列。一个阵列为X，[n_samples,n_features]作为训练样本。另一个y（注意是X和y）作为类标签，size[n_samples].
from sklearn import svm
x=[[0,0],[1,1]]
y=[0,1]
clf=svm.SVC()
#这样clf就表示SVC训练算法
clf.fit(X,y) #根据特征，进行训练，提取特征




进行训练完毕之后，可以进行预测了。


clf.predict([[2., 2.]])
免责声明：本站部分内容、图片、文字、视频等来自于互联网，仅供大家学习与交流。相关内容如涉嫌侵犯您的知识产权或其他合法权益，请向本站发送有效通知，我们会及时处理。反馈邮箱&&&&。
学生服务号
在线咨询，奖学金返现，名师点评，等你来互动SVM 支持向量机
原理就不赘述了，其余的文章有讲过。SVM是一种十分优秀的分类算法，使用SVM也能给股票进行一定程度上的预测。
因为是分类算法，因此不像ARIMA一样预测的是时序。分类就要有东西可分，因此将当日涨记为1，跌记为0，作为分类的依据。使用历史数据作为训练数据。
处理数据：
1.股票历史数据来源于yahoo_finance api，获取其中Open，Close，Low，High，Volume作为基础。因为除去Volume以外，其余数据都是Price，基于Price并不能很好的表达股票的特性，或者说并不太适用于SVM分类算法的特性。基于SVM算法的特性，股票并不是到达一个价格范围就有大概率涨或跌（不知道我这个表达大家能不能看懂）。
2.基于上述原因，我决定将Price转换成另一种形式的数据。例如：High-Low=全天最大价格差，Open-YesterdayOpen=当天Open价格变动，Open-YesterdayClose=开盘价格变动。（我也并不太懂经济学，仅仅是为了寻找另一种更好的方案）
3.单纯地基于历史数据是完全不够的，因此还使用了R语言和tm.plugin.sentiment包，进行语义分析，进行新闻正面负面的判定。这块不是我做的，了解的并不多。新闻并不是每天都有的，这样的话新闻数据就显得有些鸡肋，无法在分类算法中起到作用，但是我们能在多个站点中提取，或是直接将关键字定为Debt（判断大众反应）。
4.这里仅仅是进行了两个站点的新闻挖掘，然后可通过rpy2包在Python中运行R语言，或是R语言得到的数据导出成Json，Python再读取。至此，数据处理告一段落。
股票数据不能完全基于历史数据，因此需要一定数量的历史数据推出预测数据，例如这边使用了70天的数据训练，来推出后一天的股票涨跌，而不是所有的历史数据。
最后的成绩是53.74%的正确率，对于一个基本使用历史数据来预测股市的方法而言已经是个不错的结局了。
/jerry81333/StockProdiction/
本文已收录于以下专栏：
相关文章推荐
全称为自回归积分滑动平均模型(Autoregressive Integrated Moving Average Model,简记ARIMA)。核心函数是ARIMA（p，d，q）称为差分自回...
本文由北邮@爱可可-爱生活 老师推荐，阿里云云栖社区组织翻译。
以下为译文
本篇文章是&Python股市数据分析&两部曲中的第一部分（第二部分的文章在这里），内容基于我在犹他州立大学MAT...
刚刚开始进入量化投资领域，最近在做金融数据方面的预测，用到了数据挖掘的
摘要: 前几天列出了100个深度学习的源代码项目，但其中大部分都不活跃。这里我们精选出18个最活跃的项目，每个都制作了信息卡片，一目了然，方便比较和转贴。 ...
...
什么是机器学习 （Machine Learning）
      机器学习是研究计算机怎样模拟或实现人类的学习行为，以获取新的知识或技能，重新组织已有的知识结构使之不断改善自身的性能。它是人工智能的...
机器学习算法与Python实践之（二）支持向量机（SVM）初级http://blog.csdn.net/zouxy09        机器学习算法与Python实践这个系列...
本策略是为了验证SVM对于大盘涨跌的预测是否有效，相比于纯随机策略，是否有明显的提高。
SVM模型用06~14年的数据训练，16年1月~12月的数据用来回测，这样是为了避免因为在模型中投入了现阶段的...
SVM 支持向量机
他的最新文章
讲师：王哲涵
讲师：韦玮
您举报文章：
举报原因：
原文地址：
原因补充：
(最多只允许输入30个字)机器学习算法与Python实践之（四）支持向量机（SVM）实现(源码)
& & & &机器学习算法与Python实践这个系列主要是参考这本书。因为自己想学习Python，然后也想对一些机器学习算法加深下了解，所以就想通过Python来实现几个比较常用的机器学习算法。恰好遇见这本同样定位的书籍，所以就参考这本书的过程来学习了。
& & & &在这一节我们主要是对支持向量机进行系统的回顾，以及通过Python来实现。由于内容很多，所以这里分成三篇博文。第一篇讲SVM初级，第二篇讲进阶，主要是把SVM整条知识链理直，第三篇介绍Python的实现。SVM有很多介绍的非常好的博文，具体可以参考本文列出的参考文献和推荐阅读资料。在本文中，定位在于把集大成于一身的SVM的整体知识链理直，所以不会涉及细节的推导。网上的解说的很好的推导和书籍很多，大家可以进一步参考。
二、线性可分SVM与硬间隔最大化
三、Dual优化问题
& & & &3.1、对偶问题
& & & &3.2、SVM优化的对偶问题
四、松弛向量与软间隔最大化
五、核函数
六、多类分类之SVM
& & & &6.1、“一对多”的方法
& & & &6.2、“一对一”的方法
七、KKT条件分析
八、SVM的实现之SMO算法
& & & &8.1、坐标下降算法
& & & &8.2、SMO算法原理
& & & &8.3、SMO算法的Python实现
九、参考文献与推荐阅读
八、SVM的实现之SMO算法
& & & 终于到SVM的实现部分了。那么神奇和有效的东西还得回归到实现才可以展示其强大的功力。SVM有效而且存在很高效的训练算法，这也是工业界非常青睐SVM的原因。
& & & 前面讲到，SVM的学习问题可以转化为下面的对偶问题：
& & & &需要满足的KKT条件：
& & & &也就是说找到一组αi可以满足上面的这些条件的就是该目标的一个最优解。所以我们的优化目标是找到一组最优的αi*。一旦求出这些αi*，就很容易计算出权重向量w*和b，并得到分隔超平面了。
& & & &这是个凸二次规划问题，它具有全局最优解，一般可以通过现有的工具来优化。但当训练样本非常多的时候，这些优化算法往往非常耗时低效，以致无法使用。从SVM提出到现在，也出现了很多优化训练的方法。其中，非常出名的一个是1982年由Microsoft Research的John C. Platt在论文《Sequential Minimal Optimization: A Fast Algorithm for TrainingSupport Vector
Machines》中提出的Sequential Minimal Optimization序列最小化优化算法，简称SMO算法。SMO算法的思想很简单，它将大优化的问题分解成多个小优化的问题。这些小问题往往比较容易求解，并且对他们进行顺序求解的结果与将他们作为整体来求解的结果完全一致。在结果完全一致的同时，SMO的求解时间短很多。在深入SMO算法之前，我们先来了解下坐标下降这个算法，SMO其实基于这种简单的思想的。
8.1、坐标下降（上升）法
& & & 假设要求解下面的优化问题：
& & & 在这里，我们需要求解m个变量αi，一般来说是通过梯度下降（这里是求最大值，所以应该叫上升）等算法每一次迭代对所有m个变量αi也就是α向量进行一次性优化。通过误差每次迭代调整α向量中每个元素的值。而坐标上升法（坐标上升与坐标下降可以看做是一对，坐标上升是用来求解max最优化问题，坐标下降用于求min最优化问题）的思想是每次迭代只调整一个变量αi的值，其他变量的值在这次迭代中固定不变。
& & & &最里面语句的意思是固定除αi之外的所有αj(i不等于j)，这时W可看作只是关于αi的函数，那么直接对αi求导优化即可。这里我们进行最大化求导的顺序i是从1到m，可以通过更改优化顺序来使W能够更快地增加并收敛。如果W在内循环中能够很快地达到最优，那么坐标上升法会是一个很高效的求极值方法。
& & & 用个二维的例子来说明下坐标下降法：我们需要寻找f(x,y)=x2+xy+y2的最小值处的(x*, y*)，也就是下图的F*点的地方。
& & & &假设我们初始的点是A（图是函数投影到xoy平面的等高线图，颜色越深值越小），我们需要达到F*的地方。那最快的方法就是图中黄色线的路径，一次性就到达了，其实这个是牛顿优化法，但如果是高维的话，这个方法就不太高效了（因为需要求解矩阵的逆，这个不在这里讨论）。我们也可以按照红色所指示的路径来走。从A开始，先固定x，沿着y轴往让f(x, y)值减小的方向走到B点，然后固定y，沿着x轴往让f(x, y)值减小的方向走到C点，不断循环，直到到达F*。反正每次只要我们都往让f(x,
y)值小的地方走就行了，这样脚踏实地，一步步走，每一步都使f(x, y)慢慢变小，总有一天，皇天不负有心人的。到达F*也是时间问题。到这里你可能会说，这红色线比黄色线贫富差距也太严重了吧。因为这里是二维的简单的情况嘛。如果是高维的情况，而且目标函数很复杂的话，再加上样本集很多，那么在梯度下降中，目标函数对所有αi求梯度或者在牛顿法中对矩阵求逆，都是很耗时的。这时候，如果W只对单个αi优化很快的时候，坐标下降法可能会更加高效。
8.2、SMO算法
& & & &SMO算法的思想和坐标下降法的思想差不多。唯一不同的是，SMO是一次迭代优化两个α而不是一个。为什么要优化两个呢？
& & & &我们回到这个优化问题。我们可以看到这个优化问题存在着一个约束，也就是
& & & &假设我们首先固定除α1以外的所有参数，然后在α1上求极值。但需要注意的是，因为如果固定α1以外的所有参数，由上面这个约束条件可以知道，α1将不再是变量（可以由其他值推出），因为问题中规定了：
& & & 因此，我们需要一次选取两个参数做优化，比如αi和αj，此时αi可以由αj和其他参数表示出来。这样回代入W中，W就只是关于αj的函数了，这时候就可以只对αj进行优化了。在这里就是对αj进行求导，令导数为0就可以解出这个时候最优的αj了。然后也可以得到αi。这就是一次的迭代过程，一次迭代只调整两个拉格朗日乘子αi和αj。SMO之所以高效就是因为在固定其他参数后，对一个参数优化过程很高效（对一个参数的优化可以通过解析求解，而不是迭代。虽然对一个参数的一次最小优化不可能保证其结果就是所优化的拉格朗日乘子的最终结果，但会使目标函数向极小值迈进一步，这样对所有的乘子做最小优化，直到所有满足KKT条件时，目标函数达到最小）。
& & & &总结下来是：
重复下面过程直到收敛{
（1）选择两个拉格朗日乘子αi和αj；
（2）固定其他拉格朗日乘子αk(k不等于i和j)，只对αi和αj优化w(α);
（3）根据优化后的αi和αj，更新截距b的值；
& & & & 那训练里面这两三步骤到底是怎么实现的，需要考虑什么呢？下面我们来具体分析下：
（1）选择αi和αj：
& & & & 我们现在是每次迭代都优化目标函数的两个拉格朗日乘子αi和αj，然后其他的拉格朗日乘子保持固定。如果有N个训练样本，我们就有N个拉格朗日乘子需要优化，但每次我们只挑两个进行优化，我们就有N(N-1)种选择。那到底我们要选择哪对αi和αj呢？选择哪对才好呢？想想我们的目标是什么？我们希望把所有违法KKT条件的样本都纠正回来，因为如果所有样本都满足KKT条件的话，我们的优化就完成了。那就很直观了，哪个害群之马最严重，我们得先对他进行思想教育，让他尽早回归正途。OK，那我们选择的第一个变量αi就选违法KKT条件最严重的那一个。那第二个变量αj怎么选呢？
& & & &我们是希望快点找到最优的N个拉格朗日乘子，使得代价函数最大，换句话说，要最快的找到代价函数最大值的地方对应的N个拉格朗日乘子。这样我们的训练时间才会短。就像你从广州去北京，有飞机和绿皮车给你选，你选啥？（就算你不考虑速度，也得考虑下空姐的感受嘛，别辜负了她们渴望看到你的期盼，哈哈）。有点离题了，anyway，每次迭代中，哪对αi和αj可以让我更快的达到代价函数值最大的地方，我们就选他们。或者说，走完这一步，选这对αi和αj代价函数值增加的值最多，比选择其他所有αi和αj的结合中都多。这样我们才可以更快的接近代价函数的最大值，也就是达到优化的目标了。再例如，下图，我们要从A点走到B点，按蓝色的路线走c2方向的时候，一跨一大步，按红色的路线走c1方向的时候，只能是人类的一小步。所以，蓝色路线走两步就迈进了成功之门，而红色的路线，人生曲折，好像成功遥遥无期一样，故曰，选择比努力更重要！
& & & &真啰嗦！说了半天，其实就一句话：为什么每次迭代都要选择最好的αi和αj，就是为了更快的收敛！那实践中每次迭代到底要怎样选αi和αj呢？这有个很好听的名字叫启发式选择，主要思想是先选择最有可能需要优化（也就是违反KKT条件最严重）的αi，再针对这样的αi选择最有可能取得较大修正步长的αj。具体是以下两个过程：
1）第一个变量αi的选择：
& & & &SMO称选择第一个变量的过程为外层循环。外层训练在训练样本中选取违法KKT条件最严重的样本点。并将其对应的变量作为第一个变量。具体的，检验训练样本(xi, yi)是否满足KKT条件，也就是：
& & & &该检验是在ε范围内进行的。在检验过程中，外层循环首先遍历所有满足条件0&αj&C的样本点，即在间隔边界上的支持向量点，检验他们是否满足KKT条件，然后选择违反KKT条件最严重的αi。如果这些样本点都满足KKT条件，那么遍历整个训练集，检验他们是否满足KKT条件，然后选择违反KKT条件最严重的αi。
& & & &优先选择遍历非边界数据样本，因为非边界数据样本更有可能需要调整，边界数据样本常常不能得到进一步调整而留在边界上。由于大部分数据样本都很明显不可能是支持向量，因此对应的α乘子一旦取得零值就无需再调整。遍历非边界数据样本并选出他们当中违反KKT 条件为止。当某一次遍历发现没有非边界数据样本得到调整时，遍历所有数据样本，以检验是否整个集合都满足KKT条件。如果整个集合的检验中又有数据样本被进一步进化，则有必要再遍历非边界数据样本。这样，不停地在遍历所有数据样本和遍历非边界数据样本之间切换，直到整个样本集合都满足KKT条件为止。以上用KKT条件对数据样本所做的检验都以达到一定精度ε就可以停止为条件。如果要求十分精确的输出算法，则往往不能很快收敛。
& & & &对整个数据集的遍历扫描相当容易，而实现对非边界αi的扫描时，首先需要将所有非边界样本的αi值（也就是满足0&αi&C）保存到新的一个列表中，然后再对其进行遍历。同时，该步骤跳过那些已知的不会改变的αi值。
2）第二个变量αj的选择：
& & & &在选择第一个αi后，算法会通过一个内循环来选择第二个αj值。因为第二个乘子的迭代步长大致正比于|Ei-Ej|，所以我们需要选择能够最大化|Ei-Ej|的第二个乘子（选择最大化迭代步长的第二个乘子）。在这里，为了节省计算时间，我们建立一个全局的缓存用于保存所有样本的误差值，而不用每次选择的时候就重新计算。我们从中选择使得步长最大或者|Ei-Ej|最大的αj。
（2）优化αi和αj：
& & & &选择这两个拉格朗日乘子后，我们需要先计算这些参数的约束值。然后再求解这个约束最大化问题。
& & & &首先，我们需要给αj找到边界L&=αj&=H，以保证αj满足0&=αj&=C的约束。这意味着αj必须落入这个盒子中。由于只有两个变量(αi, αj)，约束可以用二维空间中的图形来表示，如下图：
& & & &不等式约束使得(αi,αj)在盒子[0, C]x[0, C]内，等式约束使得(αi, αj)在平行于盒子[0, C]x[0, C]的对角线的直线上。因此要求的是目标函数在一条平行于对角线的线段上的最优值。这使得两个变量的最优化问题成为实质的单变量的最优化问题。由图可以得到，αj的上下界可以通过下面的方法得到：
& & & &我们优化的时候，αj必须要满足上面这个约束。也就是说上面是αj的可行域。然后我们开始寻找αj，使得目标函数最大化。通过推导得到αj的更新公式如下：
& & & &这里Ek可以看做对第k个样本，SVM的输出与期待输出，也就是样本标签的误差。
& & & &而η实际上是度量两个样本i和j的相似性的。在计算η的时候，我们需要使用核函数，那么就可以用核函数来取代上面的内积。
& & & &得到新的αj后，我们需要保证它处于边界内。换句话说，如果这个优化后的值跑出了边界L和H，我们就需要简单的裁剪，将αj收回这个范围：
& & & &最后，得到优化的αj后，我们需要用它来计算αi：
& & & &到这里，αi和αj的优化就完成了。
（3）计算阈值b：
& & & &优化αi和αj后，我们就可以更新阈值b，使得对两个样本i和j都满足KKT条件。如果优化后αi不在边界上（也就是满足0&αi&C，这时候根据KKT条件，可以得到yigi(xi)=1，这样我们才可以计算b），那下面的阈值b1是有效的，因为当输入xi时它迫使SVM输出yi。
& & & &同样，如果0&αj&C，那么下面的b2也是有效的：
& & & 如果0&αi&C和0&αj&C都满足，那么b1和b2都有效，而且他们是相等的。如果他们两个都处于边界上（也就是αi=0或者αi=C，同时αj=0或者αj=C），那么在b1和b2之间的阈值都满足KKT条件，一般我们取他们的平均值b=(b1+b2)/2。所以，总的来说对b的更新如下：
& & & &每做完一次最小优化，必须更新每个数据样本的误差，以便用修正过的分类面对其他数据样本再做检验，在选择第二个配对优化数据样本时用来估计步长。
（4）凸优化问题终止条件：
& & & &SMO算法的基本思路是：如果说有变量的解都满足此最优化问题的KKT条件，那么这个最优化问题的解就得到了。因为KKT条件是该最优化问题的充分必要条件（证明请参考文献）。所以我们可以监视原问题的KKT条件，所以所有的样本都满足KKT条件，那么就表示迭代结束了。但是由于KKT条件本身是比较苛刻的，所以也需要设定一个容忍值，即所有样本在容忍值范围内满足KKT条件则认为训练可以结束；当然了，对于对偶问题的凸优化还有其他终止条件，可以参考文献。
8.3、SMO算法的Python实现
8.3.1、Python的准备工作
& & & 我使用的Python是2.7.5版本的。附加的库有Numpy和Matplotlib。而Matplotlib又依赖dateutil和pyparsing两个库，所以我们需要安装以上三个库。前面两个库还好安装，直接在官网下对应版本就行。但我找后两个库的时候，就没那么容易了。后来发现，其实对Python的库的下载和安装可以借助pip工具的。这个是安装和管理Python包的工具。感觉它有点像ubuntu的apt-get，需要安装什么库，直接下载和安装一条龙服务。
& & & &首先，我们需要到pip的官网：下载对应我们python版本的pip，例如我的是pip-1.4.1.tar.gz。但安装pip需要另一个工具，也就是setuptools，我们到下载ez_setup.py这个文件回来。然后在CMD命令行中执行：（注意他们的路径）
#python ez_setup.py
这时候，就会自动下载.egg等等文件然后安装完成。
& & & 然后我们解压pip-1.4.1.tar.gz。进入到该目录中，执行：
#python setup.py install
这时候就会自动安装pip到你python目录下的Scripts文件夹中。我的是C:\Python27\Scripts。
& & & &在里面我们可以看到pip.exe，然后我们进入到该文件夹中：
#cd C:\Python27\Scripts
#pip install dateutil
#pip install pyparsing
这样就可以把这些额外的库给下载回来了。非常高端大气上档次！
8.3.2、SMO算法的Python实现
& & & &在代码中已经有了比较详细的注释了。不知道有没有错误的地方，如果有，还望大家指正（每次的运行结果都有可能不同，另外，感觉有些结果似乎不太正确，但我还没发现哪里出错了，如果大家找到有错误的地方，还望大家指点下，衷心感谢）。里面我写了个可视化结果的函数，但只能在二维的数据上面使用。直接贴代码：
&&&&&&&&&&&&&&&&from&numpy&import&*&&import&time&&import&matplotlib.pyplot&as&plt&&&&&&&&&def&calcKernelValue(matrix_x,&sample_x,&kernelOption):&&&&&&kernelType&=&kernelOption[0]&&&&&&numSamples&=&matrix_x.shape[0]&&&&&&kernelValue&=&mat(zeros((numSamples,&1)))&&&&&&&&&&&&if&kernelType&==&'linear':&&&&&&&&&&kernelValue&=&matrix_x&*&sample_x.T&&&&&&elif&kernelType&==&'rbf':&&&&&&&&&&sigma&=&kernelOption[1]&&&&&&&&&&if&sigma&==&0:&&&&&&&&&&&&&&sigma&=&1.0&&&&&&&&&&for&i&in&xrange(numSamples):&&&&&&&&&&&&&&diff&=&matrix_x[i,&:]&-&sample_x&&&&&&&&&&&&&&kernelValue[i]&=&exp(diff&*&diff.T&/&(-2.0&*&sigma**2))&&&&&&else:&&&&&&&&&&raise&NameError('Not&support&kernel&type!&You&can&use&linear&or&rbf!')&&&&&&return&kernelValue&&&&&&&&def&calcKernelMatrix(train_x,&kernelOption):&&&&&&numSamples&=&train_x.shape[0]&&&&&&kernelMatrix&=&mat(zeros((numSamples,&numSamples)))&&&&&&for&i&in&xrange(numSamples):&&&&&&&&&&kernelMatrix[:,&i]&=&calcKernelValue(train_x,&train_x[i,&:],&kernelOption)&&&&&&return&kernelMatrix&&&&&&&&class&SVMStruct:&&&&&&def&__init__(self,&dataSet,&labels,&C,&toler,&kernelOption):&&&&&&&&&&self.train_x&=&dataSet&&&&&&&&&&&self.train_y&=&labels&&&&&&&&&&&&self.C&=&C&&&&&&&&&&&&&&&&&&&&&&&self.toler&=&toler&&&&&&&&&&&&&&&self.numSamples&=&dataSet.shape[0]&&&&&&&&&&&self.alphas&=&mat(zeros((self.numSamples,&1)))&&&&&&&&&&&self.b&=&0&&&&&&&&&&self.errorCache&=&mat(zeros((self.numSamples,&2)))&&&&&&&&&&self.kernelOpt&=&kernelOption&&&&&&&&&&self.kernelMat&=&calcKernelMatrix(self.train_x,&self.kernelOpt)&&&&&&&&&&&&&&&&def&calcError(svm,&alpha_k):&&&&&&output_k&=&float(multiply(svm.alphas,&svm.train_y).T&*&svm.kernelMat[:,&alpha_k]&+&svm.b)&&&&&&error_k&=&output_k&-&float(svm.train_y[alpha_k])&&&&&&return&error_k&&&&&&&&def&updateError(svm,&alpha_k):&&&&&&error&=&calcError(svm,&alpha_k)&&&&&&svm.errorCache[alpha_k]&=&[1,&error]&&&&&&&&def&selectAlpha_j(svm,&alpha_i,&error_i):&&&&&&svm.errorCache[alpha_i]&=&[1,&error_i]&&&&&&&candidateAlphaList&=&nonzero(svm.errorCache[:,&0].A)[0]&&&&&&&maxStep&=&0;&alpha_j&=&0;&error_j&=&0&&&&&&&&&&&&&&if&len(candidateAlphaList)&&&1:&&&&&&&&&&for&alpha_k&in&candidateAlphaList:&&&&&&&&&&&&&&if&alpha_k&==&alpha_i:&&&&&&&&&&&&&&&&&&&continue&&&&&&&&&&&&&&error_k&=&calcError(svm,&alpha_k)&&&&&&&&&&&&&&if&abs(error_k&-&error_i)&&&maxStep:&&&&&&&&&&&&&&&&&&maxStep&=&abs(error_k&-&error_i)&&&&&&&&&&&&&&&&&&alpha_j&=&alpha_k&&&&&&&&&&&&&&&&&&error_j&=&error_k&&&&&&&&&&&&else:&&&&&&&&&&&&&&&&&&&&&alpha_j&=&alpha_i&&&&&&&&&&while&alpha_j&==&alpha_i:&&&&&&&&&&&&&&alpha_j&=&int(random.uniform(0,&svm.numSamples))&&&&&&&&&&error_j&=&calcError(svm,&alpha_j)&&&&&&&&&&&&return&alpha_j,&error_j&&&&&&&&def&innerLoop(svm,&alpha_i):&&&&&&error_i&=&calcError(svm,&alpha_i)&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&if&(svm.train_y[alpha_i]&*&error_i&&&-svm.toler)&and&(svm.alphas[alpha_i]&&&svm.C)&or\&&&&&&&&&&(svm.train_y[alpha_i]&*&error_i&&&svm.toler)&and&(svm.alphas[alpha_i]&&&0):&&&&&&&&&&&&&&&&&&&&&&alpha_j,&error_j&=&selectAlpha_j(svm,&alpha_i,&error_i)&&&&&&&&&&alpha_i_old&=&svm.alphas[alpha_i].copy()&&&&&&&&&&alpha_j_old&=&svm.alphas[alpha_j].copy()&&&&&&&&&&&&&&&&&&&&&&if&svm.train_y[alpha_i]&!=&svm.train_y[alpha_j]:&&&&&&&&&&&&&&L&=&max(0,&svm.alphas[alpha_j]&-&svm.alphas[alpha_i])&&&&&&&&&&&&&&H&=&min(svm.C,&svm.C&+&svm.alphas[alpha_j]&-&svm.alphas[alpha_i])&&&&&&&&&&else:&&&&&&&&&&&&&&L&=&max(0,&svm.alphas[alpha_j]&+&svm.alphas[alpha_i]&-&svm.C)&&&&&&&&&&&&&&H&=&min(svm.C,&svm.alphas[alpha_j]&+&svm.alphas[alpha_i])&&&&&&&&&&if&L&==&H:&&&&&&&&&&&&&&return&0&&&&&&&&&&&&&&&&&&&&&&eta&=&2.0&*&svm.kernelMat[alpha_i,&alpha_j]&-&svm.kernelMat[alpha_i,&alpha_i]&\&&&&&&&&&&&&&&&&&&&&-&svm.kernelMat[alpha_j,&alpha_j]&&&&&&&&&&if&eta&&=&0:&&&&&&&&&&&&&&return&0&&&&&&&&&&&&&&&&&&&&&&svm.alphas[alpha_j]&-=&svm.train_y[alpha_j]&*&(error_i&-&error_j)&/&eta&&&&&&&&&&&&&&&&&&&&&&if&svm.alphas[alpha_j]&&&H:&&&&&&&&&&&&&&svm.alphas[alpha_j]&=&H&&&&&&&&&&if&svm.alphas[alpha_j]&&&L:&&&&&&&&&&&&&&svm.alphas[alpha_j]&=&L&&&&&&&&&&&&&&&&&&&&&&if&abs(alpha_j_old&-&svm.alphas[alpha_j])&&&0.00001:&&&&&&&&&&&&&&updateError(svm,&alpha_j)&&&&&&&&&&&&&&return&0&&&&&&&&&&&&&&&&&&&&&&svm.alphas[alpha_i]&+=&svm.train_y[alpha_i]&*&svm.train_y[alpha_j]&\&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*&(alpha_j_old&-&svm.alphas[alpha_j])&&&&&&&&&&&&&&&&&&&&&&b1&=&svm.b&-&error_i&-&svm.train_y[alpha_i]&*&(svm.alphas[alpha_i]&-&alpha_i_old)&\&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*&svm.kernelMat[alpha_i,&alpha_i]&\&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&-&svm.train_y[alpha_j]&*&(svm.alphas[alpha_j]&-&alpha_j_old)&\&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*&svm.kernelMat[alpha_i,&alpha_j]&&&&&&&&&&b2&=&svm.b&-&error_j&-&svm.train_y[alpha_i]&*&(svm.alphas[alpha_i]&-&alpha_i_old)&\&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*&svm.kernelMat[alpha_i,&alpha_j]&\&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&-&svm.train_y[alpha_j]&*&(svm.alphas[alpha_j]&-&alpha_j_old)&\&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*&svm.kernelMat[alpha_j,&alpha_j]&&&&&&&&&&if&(0&&&svm.alphas[alpha_i])&and&(svm.alphas[alpha_i]&&&svm.C):&&&&&&&&&&&&&&svm.b&=&b1&&&&&&&&&&elif&(0&&&svm.alphas[alpha_j])&and&(svm.alphas[alpha_j]&&&svm.C):&&&&&&&&&&&&&&svm.b&=&b2&&&&&&&&&&else:&&&&&&&&&&&&&&svm.b&=&(b1&+&b2)&/&2.0&&&&&&&&&&&&&&&&&&&&&&updateError(svm,&alpha_j)&&&&&&&&&&updateError(svm,&alpha_i)&&&&&&&&&&&&return&1&&&&&&else:&&&&&&&&&&return&0&&&&&&&&def&trainSVM(train_x,&train_y,&C,&toler,&maxIter,&kernelOption&=&('rbf',&1.0)):&&&&&&&&&&&&startTime&=&time.time()&&&&&&&&&&&&&&svm&=&SVMStruct(mat(train_x),&mat(train_y),&C,&toler,&kernelOption)&&&&&&&&&&&&&&&&&&entireSet&=&True&&&&&&alphaPairsChanged&=&0&&&&&&iterCount&=&0&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&while&(iterCount&&&maxIter)&and&((alphaPairsChanged&&&0)&or&entireSet):&&&&&&&&&&alphaPairsChanged&=&0&&&&&&&&&&&&&&&&&&&&&&if&entireSet:&&&&&&&&&&&&&&for&i&in&xrange(svm.numSamples):&&&&&&&&&&&&&&&&&&alphaPairsChanged&+=&innerLoop(svm,&i)&&&&&&&&&&&&&&print&'---iter:%d&entire&set,&alpha&pairs&changed:%d'&%&(iterCount,&alphaPairsChanged)&&&&&&&&&&&&&&iterCount&+=&1&&&&&&&&&&&&&&&&&&&&else:&&&&&&&&&&&&&&nonBoundAlphasList&=&nonzero((svm.alphas.A&&&0)&*&(svm.alphas.A&&&svm.C))[0]&&&&&&&&&&&&&&for&i&in&nonBoundAlphasList:&&&&&&&&&&&&&&&&&&alphaPairsChanged&+=&innerLoop(svm,&i)&&&&&&&&&&&&&&print&'---iter:%d&non&boundary,&alpha&pairs&changed:%d'&%&(iterCount,&alphaPairsChanged)&&&&&&&&&&&&&&iterCount&+=&1&&&&&&&&&&&&&&&&&&&&&&if&entireSet:&&&&&&&&&&&&&&entireSet&=&False&&&&&&&&&&elif&alphaPairsChanged&==&0:&&&&&&&&&&&&&&entireSet&=&True&&&&&&&&print&'Congratulations,&training&complete!&Took&%fs!'&%&(time.time()&-&startTime)&&&&&&return&svm&&&&&&&&def&testSVM(svm,&test_x,&test_y):&&&&&&test_x&=&mat(test_x)&&&&&&test_y&=&mat(test_y)&&&&&&numTestSamples&=&test_x.shape[0]&&&&&&supportVectorsIndex&=&nonzero(svm.alphas.A&&&0)[0]&&&&&&supportVectors&&&&&&=&svm.train_x[supportVectorsIndex]&&&&&&supportVectorLabels&=&svm.train_y[supportVectorsIndex]&&&&&&supportVectorAlphas&=&svm.alphas[supportVectorsIndex]&&&&&&matchCount&=&0&&&&&&for&i&in&xrange(numTestSamples):&&&&&&&&&&kernelValue&=&calcKernelValue(supportVectors,&test_x[i,&:],&svm.kernelOpt)&&&&&&&&&&predict&=&kernelValue.T&*&multiply(supportVectorLabels,&supportVectorAlphas)&+&svm.b&&&&&&&&&&if&sign(predict)&==&sign(test_y[i]):&&&&&&&&&&&&&&matchCount&+=&1&&&&&&accuracy&=&float(matchCount)&/&numTestSamples&&&&&&return&accuracy&&&&&&&&def&showSVM(svm):&&&&&&if&svm.train_x.shape[1]&!=&2:&&&&&&&&&&print&&Sorry!&I&can&not&draw&because&the&dimension&of&your&data&is&not&2!&&&&&&&&&&&return&1&&&&&&&&&&&&&&for&i&in&xrange(svm.numSamples):&&&&&&&&&&if&svm.train_y[i]&==&-1:&&&&&&&&&&&&&&plt.plot(svm.train_x[i,&0],&svm.train_x[i,&1],&'or')&&&&&&&&&&elif&svm.train_y[i]&==&1:&&&&&&&&&&&&&&plt.plot(svm.train_x[i,&0],&svm.train_x[i,&1],&'ob')&&&&&&&&&&&&&&supportVectorsIndex&=&nonzero(svm.alphas.A&&&0)[0]&&&&&&for&i&in&supportVectorsIndex:&&&&&&&&&&plt.plot(svm.train_x[i,&0],&svm.train_x[i,&1],&'oy')&&&&&&&&&&&&&&&&&&w&=&zeros((2,&1))&&&&&&for&i&in&supportVectorsIndex:&&&&&&&&&&w&+=&multiply(svm.alphas[i]&*&svm.train_y[i],&svm.train_x[i,&:].T)&&&&&&&min_x&=&min(svm.train_x[:,&0])[0,&0]&&&&&&max_x&=&max(svm.train_x[:,&0])[0,&0]&&&&&&y_min_x&=&float(-svm.b&-&w[0]&*&min_x)&/&w[1]&&&&&&y_max_x&=&float(-svm.b&-&w[0]&*&max_x)&/&w[1]&&&&&&plt.plot([min_x,&max_x],&[y_min_x,&y_max_x],&'-g')&&&&&&plt.show()&&
& & & &测试的数据来自。有100个样本，每个样本两维，最后是对应的标签，例如：
3..977398 &&&&&&&& -1
3..556416 &&&&&&&& -1
7..580030&&&&&&&& 1
2..004466&&&&&&&& -1
& & & &测试代码中首先加载这个数据库，然后用前面80个样本来训练，再用剩下的20个样本的测试，并显示训练后的模型和分类结果。测试代码如下：
test_SVM.py
&&&&&&&&&&&&&&&&from&numpy&import&*&&import&SVM&&&&&&&&print&&step&1:&load&data...&&&dataSet&=&[]&&labels&=&[]&&fileIn&=&open('E:/Python/Machine&Learning&in&Action/testSet.txt')&&for&line&in&fileIn.readlines():&&&&&&lineArr&=&line.strip().split('\t')&&&&&&dataSet.append([float(lineArr[0]),&float(lineArr[1])])&&&&&&labels.append(float(lineArr[2]))&&&&dataSet&=&mat(dataSet)&&labels&=&mat(labels).T&&train_x&=&dataSet[0:81,&:]&&train_y&=&labels[0:81,&:]&&test_x&=&dataSet[80:101,&:]&&test_y&=&labels[80:101,&:]&&&&&&print&&step&2:&training...&&&C&=&0.6&&toler&=&0.001&&maxIter&=&50&&svmClassifier&=&SVM.trainSVM(train_x,&train_y,&C,&toler,&maxIter,&kernelOption&=&('linear',&0))&&&&&&print&&step&3:&testing...&&&accuracy&=&SVM.testSVM(svmClassifier,&test_x,&test_y)&&&&&&print&&step&4:&show&the&result...&&&&&print&'The&classify&accuracy&is:&%.3f%%'&%&(accuracy&*&100)&&SVM.showSVM(svmClassifier)&&
运行结果如下：
step&1:&load&data...&&step&2:&training...&&---iter:0&entire&set,&alpha&pairs&changed:8&&---iter:1&non&boundary,&alpha&pairs&changed:7&&---iter:2&non&boundary,&alpha&pairs&changed:1&&---iter:3&non&boundary,&alpha&pairs&changed:0&&---iter:4&entire&set,&alpha&pairs&changed:0&&Congratulations,&training&complete!&Took&0.058000s!&&step&3:&testing...&&step&4:&show&the&result...&&The&classify&accuracy&is:&100.000%&&
训练好的模型图：
&span style=&font-size:12&&#coding=utf-8
'''
Created on Nov 4, 2010
Chapter 5 source file for Machine Learing in Action
@author: Peter
支持向量机SVM
序列最小优化算法
基于最大间隔分割数据
优点 泛化错误率低 计算开销不大 结果易于解释
缺点 对参数调节和核函数的选择敏感 原始分类器不加修改仅适用于处理二类问题
适用数据类型 数值型和标称型数据
数据线性可分
分割超平面
我们希望找到离分隔超平面最
这里点到分隔面的距离被称为间隔①( margin )。我
们希望间隔尽可能地大，这是因为如果我们犯错或者在有限数据上训练分类器的话，我们希望分类器尽可能健壮
支持向量(support, vector)就是离分隔超平面最近的那些点。接下来要试着最大化支持向量
到分隔面的距离，需要找到此问题的优化求解方法。
分类器求解的优化问题输入数据黑分类器会输出一个类别标签，相当于一个类似于sigmoid函数在作用。下面将使用类似海维赛德阶跃函数(即单位阶跃函数)的函数对
w^T*x+b作用得到f(w^T*x+b) 其中当u&0时f(u)输出-1，反之则输出+1
SVM一般流程如下：
O)收集数据:可以使用任意方法。
(2)准备数据:需要数值型数据。
(3)分析数据:有助于可视化分隔超平面。
(4)训练算法:SVM的大部分时间都源自训练，该过程主要实现两个参数的调优。
(5)测试算法:十分简单的计算过程就可以实现。
(6)使用算法:几乎所有分类问题都可以使用SVM，值得一提的是，SVM本身是一个二类
分类器，对多类问题应用SVM需要对代码做一些修改。
SMO高效优化算法（序列最小优化）
Platt的SMO算法是将大优化问题分解为多个小优化问题来
求解的。这些小优化问题往往很容易求解，并且对它们进行顺序求解的结果与将它们作为整体来
求解的结果是完全一致的。在结果完全相同的同时，SMO算法的求解时间短很多。
SMO算法的目标是求出一系列alpha和b，一旦求出了这些alpha，就很容易计算出权重向量w
并得到分隔超平面。
SMO算法的工作原理是:每次循环中选择两个alpha进行优化处理。一旦找到一对合适的
alpha，那么就增大其中一个同时减小另一个。这里所谓的“合适”就是指两个alpha必须要符合
一定的条件，条件之一就是这两个alpha必须要在间隔边界之外，而其第二个条件则是这两个alpha
还没有进行过区间化处理或者不在边界上。
核方法或者说核技巧会将数据(有时是非线性数据)从一个低维空间映射到一个高维空间，
可以将一个在低维空间中的非线性问题转换成高维空间下的线性问题来求解。核方法不止在SVM
中适用，还可以用于其他算法中。而其中的径向基函数是一个常用的度量两个向量距离的核函数。
支持向量机是一个二类分类器。当用其解决多类问题时，则需要额外的方法对其进行扩展。
SVM的效果也对优化参数和所用核函数中的参数敏感。
'''
from numpy import *
from time import sleep
'''
由于改变一个alpha可能会导致该约束条件失效，因此我们总是同时改变两个alphao
为此，我们将构建一个辅助函数，用于在某个区间范围内随机选择一个整数。同时，我们也
需要另一个辅助函数，用于在数值太大时对其进行调整。
SMO算法中的辅助函数
'''
def loadDataSet(fileName):
dataMat = []; labelMat = []
fr = open(fileName)
for line in fr.readlines():
lineArr = line.strip().split('\t')
dataMat.append([float(lineArr[0]), float(lineArr[1])])
labelMat.append(float(lineArr[2]))
return dataMat,labelMat
'''
下一个函数。electJrand()有两个参数值，其中i是第一个alpha的下标，m是所有alpha的数
目。只要函数值不等于输人值i，函数就会进行随机选择。
最后一个辅助函数就是clipAlpha()，它是用于调整大于H或小于L的a扣ha值。尽管上述3
个辅助函数本身做的事情不多，但在分类器中却很有用处。
可以看得出来，这里采用的类别标签是一1和1，而不是0和to
'''
def selectJrand(i,m):
#i第一个alpha的下标
m所有alpha的数目，函数值不等于输入值i，函数随机选择
j=i #we want to select any J not equal to i
while (j==i):
j = int(random.uniform(0,m))
def clipAlpha(aj,H,L):
if aj & H:
if L & aj:
'''
该SMO函数的伪代码大致如下:
创建一个alpha向量并将其初始化为0向量
当迭代次数小于最大迭代次数时(外循环)
对数据集中的每个数据向量(内循环):
如果该数据向量可以被优化:
随机选择另外一个数据向量
同时优化这两个向量
如果两个向量都不能被优化，退出内循环
如果所有向量都没被优化，增加迭代数目，继续下一次循环
在Python中，如果某行以\符号结束，那
么就意味着该行语句没有结束并会在下一行延续。下面的代码当中有很多很长的语句必须要分成
多行来写。因此，下面的程序中使用了多个\符号。
这个函数比较大，或许是我所知道的本书中最大的一个函数。该函数有5个输人参数，分别
是:数据集、类别标签、常数C、容错率和取消前最大的循环次数。在本书，我们构建函数时采
用了通用的接口，这样就可以对算法和数据源进行组合或配对处理。上述函数将多个列表和输人
参数转换成NumPy矩阵，这样就可以简化很多数学处理操作。由于转置了类别标签，因此我们得
到的就是一个列向量而不是列表。于是类别标签向量的每行元素都和数据矩阵中的行一一对应。
我们也可以通过矩阵dataMatIn的shape属性得到常数m和n。最后，我们就可以构建一个alpha列
矩阵，矩阵中元素都初始化为。，并建立一个iter变量。该变量存储的则是在没有任何alpha改变
的情况下遍历数据集的次数。当该变量达到输人值max工ter时，函数结束运行并退出。
每次循环当中，将alphaPairsChanged先设为。，然后再对整个集合顺序遍历。变量
alphaPairsChanged用于记录alpha是否已经进行优化。当然，在循环结束时就会得知这一点。
首先，fXi能够计算出来，这就是我们预测的类别。然后，基于这个实例的预测结果和真实结果
的比对，就可以计算误差Ei。如果误差很大，那么可以对该数据实例所对应的alph植进行优化。
对该条件的测试处于上述程序清单的0处。在if语句中，不管是正间隔还是负间隔都会被测试。
并且在该if语句中，也要同时检查alpha值，以保证其不能等于0或C。由于后面alphai}“于0或大于
C时将被调整为。或C，所以一旦在该i晤句中它们等于这两个值的话，那么它们就已经在“边界”
上了，因而不再能够减小或增大，因此也就不值得再对它们进行优化了。’
接下来，可以利用程序清单6-1中的辅助函数来随机选择第二个alpha值，即alpha [ j ].。
同样，可以采用第一个alpha (alpha [ i ] )的误差计算方法，来计算这个alpha值的误差。这个
过程可以通过c opy ( )的方法来实现，因此稍后可以将新的alph植与老的alph植进行比较。Python
则会通过引用的方式传递所有列表，所以必须明确地告知Python要为alpha工。1d和alphaJold
分配新的内存;否则的话，在对新值和旧值进行比较时，我们就看不到新旧值的变化。之后我们
开始计算L和H)，它们用于将alpha[j]调整到。到C之间。如果L和H相等，就不做任何改变，直
接执行continue语句。这在Python中，则意味着本次循环结束直接运行下一次for的循环。
Eta是alpha [ j]的最优修改量，在那洲反长的计算代码行中得到。如果eta为。，那就是说
需要退出for循环的当前迭代过程。该过程对真实SMO算法进行了简化处理。如果eta为0，那么
计算新的alpha[j]就比较麻烦了，这里我们就不对此进行详细的介绍了。有需要的读者可以阅
读Platt的原文来了解更多的细节。现实中，这种情况并不常发生，因此忽略这一部分通常也无伤
大雅。于是，可以计算出一个新的alpha[j]，然后利用程序清单6-1中的辅助函数以及L与H值对
其进行调整。
然后，就是需要检查alpha [ j]是否有轻微改变。如果是的话，就退出for循环。然后，
alpha[i]和alpha[j]同样进行改变，虽然改变的大小一样，但是改变的方向正好相反(即如果
一个增加，那么另外一个减少)0。在对alpha[i]和alpha[j]进行优化之后，给这两个alpha
值设置一个常数项b)o
最后，在优化过程结束的同时，必须确保在合适的时机结束循环。如果程序执行到for循环
的最后一行都不执行continue语句，那么就已经成功地改变了一对alpha，同时可以增加
alphaPairsChanged的值。在for循环之外，需要检查alph植是否做了更新，如果有更新则将
iter设为0后继续运行程序。只有在所有数据集上遍历max工七er次，且不再发生任何alph谁改之
后，程序才会停止并退出while循环。
&&& alphas[alphas&0]
matrix([[ 0.,
由于SMO算法的随机性，读者运行后所得到的结果可能会与上述结果不同。
alphas[alphas&0〕命令是数组过滤(array filtering)的一个实例，而且它只又栩umPy类型有用，
却并不适用于Python中的正则表(regular list )。如果输人alpha&0，那么就会得到一个布尔数组，
并且在不等式成立的情况下，其对应值为正确的。于是，在将该布尔数组应用到原始的矩阵当中
时，就会得到一个NumPy矩阵，并且其中矩阵仅仅包含大于。的值。
支持向量机的个数：
&&& shape(alphas[alphas&0])
了解哪些数据点支持向量，输入：
&&& for i in range(100):
if alphas[i]&0.0:print dataArr[i],labelArr[i]
得到如下结果：
[4..507396] -1.0
[3.223038, -0.552392] -1.0
[3.457096, -0.082216] -1.0
[2.893743, -1.643468] -1.0
简化版SMO算法
'''
def smoSimple(dataMatIn, classLabels, C, toler, maxIter):
#五个输入函数，数据集 类别标签 常数C 容错率 取消前最大的循环次数
dataMatrix = mat(dataMatIn); labelMat = mat(classLabels).transpose()
b = 0; m,n = shape(dataMatrix)
#得到常数项n m
alphas = mat(zeros((m,1)))
#构建alpha列矩阵，矩阵中元素都初始化为0，并建立一个iter变量
while (iter & maxIter):
#如果alpha可以更改 进入优化过程
当变量iter达到输入值maxItem时，函数结束运行并退出
alphaPairsChanged = 0
#每次循环中 将alphaPairsChanged先设为0，然后再对整个结婚顺序遍历
alphaPairsChanged用于记录alpha是否已经进行优化
for i in range(m):
fXi = float(multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[i,:].T)) + b
#能够计算出来，这就是预测的类别
Ei = fXi - float(labelMat[i])#if checks if an example violates KKT conditions
实例的预测结果与真是的预测结果之间的比对
if ((labelMat[i]*Ei & -toler) and (alphas[i] & C)) or ((labelMat[i]*Ei & toler) and (alphas[i] & 0)):
#如果误差很大，那么可以对该数据实例所对应的alpha值进行优化
不管是正间隔还是负间隔都会被检测
同时检查alpha的值不能等于0 or C
j = selectJrand(i,m)
#随机选择第二个alpha
fXj = float(multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[j,:].T)) + b
Ej = fXj - float(labelMat[j])
alphaIold = alphas[i].copy(); alphaJold = alphas[j].copy();
if (labelMat[i] != labelMat[j]):
#保证alpha在0~c之间
L = max(0, alphas[j] - alphas[i])
#两个alpha值进行比较
H = min(C, C + alphas[j] - alphas[i])
L = max(0, alphas[j] + alphas[i] - C)
H = min(C, alphas[j] + alphas[i])
if L==H: print &L==H&; continue
#如果L 和 H的值相等 不做任何改变 直到continue
eta = 2.0 * dataMatrix[i,:]*dataMatrix[j,:].T - dataMatrix[i,:]*dataMatrix[i,:].T - dataMatrix[j,:]*dataMatrix[j,:].T
#eta最优修改量
eta为0，退出for循环的当前迭代过程
if eta &= 0: print &eta&=0&; continue
alphas[j] -= labelMat[j]*(Ei - Ej)/eta
alphas[j] = clipAlpha(alphas[j],H,L)
#检查是否有轻微改变，如果是的话，就退出for循环 然后alpha[i]和alpha[j]同样进行改变，改变的大小一样，方向相反
if (abs(alphas[j] - alphaJold) & 0.00001): print &j not moving enough&; continue
alphas[i] += labelMat[j]*labelMat[i]*(alphaJold - alphas[j])#update i by the same amount as j
对i进行修改，修改量与j相同但是方向相反
#the update is in the oppostie direction
b1 = b - Ei- labelMat[i]*(alphas[i]-alphaIold)*dataMatrix[i,:]*dataMatrix[i,:].T - labelMat[j]*(alphas[j]-alphaJold)*dataMatrix[i,:]*dataMatrix[j,:].T
b2 = b - Ej- labelMat[i]*(alphas[i]-alphaIold)*dataMatrix[i,:]*dataMatrix[j,:].T - labelMat[j]*(alphas[j]-alphaJold)*dataMatrix[j,:]*dataMatrix[j,:].T
#设置常数项
if (0 & alphas[i]) and (C & alphas[i]): b = b1
elif (0 & alphas[j]) and (C & alphas[j]): b = b2
else: b = (b1 + b2)/2.0
alphaPairsChanged += 1
#如果程序执行到for循环的最后一行都不执行continue语句，那么就已经成功地改变了一对alpha值，同事增加alphaPairsChanged的值
print &iter: %d i:%d, pairs changed %d& % (iter,i,alphaPairsChanged)
if (alphaPairsChanged == 0): iter += 1
#检查alpha是否做了更新，如果更新则将iter设为0后继续执行程序。遍历所有数据集maxIter次且不再发生任何alpha修改之后，停止并推出while循环
else: iter = 0
print &iteration number: %d& % iter
return b,alphas
'''
svlvl优化中一个特别好的地方就是，所有的运算都可以写成内积(inner product，也称点积:
的形式。向量的内积指的是两个向量相乘，之后得到单个标量或者数值。我们可以把内积运算替
换成核函数，而不必做简化处理。将内积替换成核函数的方式被称为核技巧(kernel trick)或者
核“变电”( kernel substation )
径向基函数:
径向基函数是SVM中常用的一个核函数。径向基函数是一个采用向量作为自变量的函数，能
够基于向量距离运算输出一个标量。这个距离可以是从&&0,0&向量或者其他向量开始计算的距离。
kTup是一个包含核函数信息的元组，待会儿我们就能看到它的作用
了。在初始化方法结束时，矩阵K先被构建，然后再通过调用函数kernelTrans()进行填充。全
局的K值只需计算一次。然后，当想要使用核函数时，就可以对它进行调用。这也省去了很多冗
余的计算开销。
当计算矩阵K时，该过程多次调用了函数kernelTrans()。该函数有3个输人参数:2个数值
型变量和1个元组。元组 kTup给出的是核函数的信息。元组的第一个参数是描述所用核函数类型
的一个字符串，其他2个参数则都是核函数可能需要的可选参数。该函数首先构建出了一个列向
量，然后检查元组以确定核函数的类型。这里只给出了2种选择，但是依然可以很容易地通过添
加elif语句来扩展到更多选项。
在线性核函数的情况下，内积计算在“所有数据集”和“数据集中的一行”这两个输人之间
展开。在径向基核函数的情况下，在for循环中对于矩阵的每个元素计算高斯函数的值。而在for
循环结束之后，我们将计算过程应用到整个向量上去。值得一提的是，在NumPy矩阵中，除法符
号意味着对矩阵元素展开计算而不像在MATLAB中一样计算矩阵的逆
最后，如果遇到一个无法识别的元组，程序就会抛出异常，因为在这种情况下不希望程序再
继续运行，这一点相当重要。
核转换函数
'''
def kernelTrans(X, A, kTup): #calc the kernel or transform data to a higher dimensional space
X,A数值型变量 kTup元祖核函数的信息 元祖第一个参数描述核函数类型的一个字符串
其他两个参数都是核函数可能需要的可选参数
m,n = shape(X)
K = mat(zeros((m,1)))
if kTup[0]=='lin': K = X * A.T
#linear kernel
elif kTup[0]=='rbf':
for j in range(m):
deltaRow = X[j,:] - A
K[j] = deltaRow*deltaRow.T
K = exp(K/(-1*kTup[1]**2)) #divide in NumPy is element-wise not matrix like Matlab 元素间作除法
else: raise NameError('Houston We Have a Problem -- \
That Kernel is not recognized')
#遇到无法是别的元祖 抛出异常
'''
利用完整Platt SMO算法加速优化
在几百个点组成的小规模数据集上，简化版SMO算法的运行是没有什么问题的，但是在更大
的数据集上的运行速度就会变慢。刚才已经讨论了简化版SMO算法，下面我们就讨论完整版的
Platt SMO算法。在这两个版本中，实现alpha的更改和代数运算的优化环节一模一样。在优化过
程中，唯一的不同就是选择alpha的方式。完整版的Platt SMO算法应用了一些能够提速的启发方fa
Platt SMO算法是通过一个外循环来选择第一个a如h植的，并且其选择过程会在两种方式之
间进行交替:一种方式是在所有数据集上进行单遍扫描，另一种方式则是在非边界alpha中实现单
遍扫描。而所谓非边界alpha指的就是那些不等于边界。或C的alpha值。对整个数据集的扫描相当
容易，而实现非边界alpha值的扫描时，首先需要建立这些alpha值的列表，然后再对这个表进行
遍历。同时，该步骤会跳过那些已知的不会改变的a扣ha值。
在选择第一个alpha值后，算法会通过一个内循环来选择第二个alpha值。在优化过程中，会
通过最大化步长的方式来获得第二个alpha值。在简化版SMO算法中，我们会在选择j之后计算错
误率Ej。但在这里，我们会建立一个全局的缓存用于保存误差值，并从中选择使得步长或者说
Ei一Ej最大的alphas 值
首要的事情就是建立一个数据结构来保存所有的重要值，而这个过程可以通过一个对象来完
成。这里使用对象的目的并不是为了面向对象的编程，而只是作为一个数据结构来使用对象。在
将值传给函数时，我们可以通过将所有数据移到一个结构中来实现，这样就可以省掉手工输人的
麻烦了。而此时，数据就可以通过一个对象来进行传递。实际上，当完成其实现时，可以很容易
通过Python的字典来完成。但是在访问对象成员变量时，这样做会有更多的手工输人操作，对比
一下myObject.X和myObject['X']就可以知道这一点。为达到这个目的，需要构建一个仅包含
ini七方法的。ptStruc七类。该方法可以实现其成员变量的填充。除了增加了一个mx2的矩阵成
员变量eCache之外.，这些做法和简化版SMO一模一样。eCache的第一列给出的是eCache是
否有效的标志位，而第二列给出的是实际的E值。
对于给定的alpha，第一个辅助函数calcEk()能够计算E值并返回。以前，该过程是采用内嵌
的方式来完成的，但是由于该过程在这个版本的SMO算法中出现频繁，这里必须要将其单别I拎出来。
SMO算法是通过一个外循环来选择第一个alpha值的，并且其选择过程会在两种方式间进行交替
一种方式是在所有数据集上进行单遍扫描 另一种方式是在非边界alpha上进行一遍扫描
所谓的非边界alpha指的是那些不等于0或者C的alpha值
对整个数据集的扫描相当容易，而实现非边界alpha扫描时
首先需要建立这些alpha值的列表，然后再对这个表进行遍历
同时 该步骤会跳过那些已知的不会改变的alpha值
在选择第一个alpha值后，算法会通过一个内循环来选择第二个alpha。在优化过程中，会通过最大化步长的
的方式获得第二个alpha值。建立一个全局的缓存，用于保存误差值
并从中选择时的步长(Ei-Ej)最大的值
'''
class optStruct:
#建立一个数据结构保存所有重要值
def __init__(self,dataMatIn, classLabels, C, toler, kTup):
# Initialize the structure with the parameters
kTup包含核函数的元祖
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
self.K = mat(zeros((self.m,self.m)))
#构建空矩阵K 之后调用kernelTrans函数进行填充
全局K只需要计算一次
for i in range(self.m):
self.K[:,i] = kernelTrans(self.X, self.X[i,:], kTup)
def calcEk(oS, k):
#计算E值并返回
fXk = float(multiply(oS.alphas,oS.labelMat).T*oS.K[:,k] + oS.b)
Ek = fXk - float(oS.labelMat[k])
'''
下一个函数selec七J()用于选择第二个alpha或者说内循环的alpha值4。回想一下，这里的
目标是选择合适的第二个alpha值以保证在每次优化中采用最大步长。该函数的误差值与第一个
alpha}Ei和下标i有关。首先将输人值Ei在缓存中设置成为有效的。这里的有效(valid)意味着
它已经计算好了。在eCache中，代码nonzero(oS.eCache[:,0].A)[0]构建出了一个非零表。
NumPy函数nonzero()返回了一个列表，而这个列表中包含以输人列表为目录的列表值，当然
读者可以猜得到，这里的值并非零。nonzero()语句返回的是非零E值所对应的alpha值，而不是
E值本身。程序会在所有的值上进行循环并选择其中使得改变最大的那个值O。如果这是第一次
循环的话，那么就随机选择一个alpha值。当然，也存在有许多更复杂的方式来处理第一次循环的
I清况，而上述做法就能够满足我们的目的。
程序清单6-3的最后一个辅助函数是updateEk()，它会计算误差值并存人缓存当中。在对
alpha进行优化之后会用到这个值。
程序清单6-3中的代码本身的作用并不大，但是当和优化过程及外循环组合在一起时，就能
组成强大的SMO算法。
'''
def selectJ(i, oS, Ei):
#this is the second choice -heurstic, and calcs Ej
内循环的启发式方法
选择第二个alpha(内循环的alpha值)
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]
#构建一个非0表
返回一个列表
包含以输入列表为目录的列表值
E值对应的alpha值
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 = maxDeltaE = deltaE; Ej = Ek
#选择具有最大步长的j值
对所有值进行循环并获得改变最大 的那个值
第一次循环随机选取alpha值
return maxK, Ej
#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
计算误差值并存入缓存中，对alpha值进行优化之后用到这个值
Ek = calcEk(oS, k)
oS.eCache[k] = [1,Ek]
'''
完整的Platt SMO算法中的优化例程
程序清单6-4中的代码几乎和程序清单6-2中给出的smoSimple()函数一模一样，但是这里的
代码已经使用了自己的数据结构。该结构在参数。S中传递。第二个重要的修改就是使用程序清单
6-3中的SelectJ()而不是selectJrand()来选择第二个alpha的值.。最后，在alpha值改变时
更新Ecache)。
'''
def innerL(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
第二个alpha选择中的启发式方法
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])
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
else: return 0
'''
完整版Platt SMO的外循环代码
程序清单6-5给出的是完整版的Platt SMO算法，其输人和函数smoSimple()完全一样。函数
一开始构建一个数据结构来容纳所有的数据，然后需要对控制函数退出的一些变量进行初始化。
整个代码的主体是while循环，这与。moSimple()有些类似，但是这里的循环退出条件更多一些。
当迭代次数超过指定的最大值，或者遍历整个集合都未对任意alpha对进行修改时，就退出循环。
这里的max工七er变量和函数smoSimple()中的作用有一点不同，后者当没有任何alpha发生改变
时会将整个集合的一次遍历过程计成一次迭代，而这里的一次迭代定义为一次循环过程，而不管
该循环具体做了什么事。此时，如果在优化过程中存在波动就会停止，因此这里的做法优于
smoSimple()函数中的计数方法。
while循环的内部与smoSimple()中有所不同，一开始的for循环在数据集上遍历任意可能
的alpha0。我们通过调用innerL()来选择第二个alpha，并在可能时对其进行优化处理。如果有
任意一对alpha值发生改变，那么会返回1。第二个for循环遍历所有的非边界alpha值，也就是不
在边界0或C上的值.。
接下来，我们对for循环在非边界循环和完整遍历之间进行切换，并打印出迭代次数。最后
程序将会返回常数b和alpha值。
'''
def smoP(dataMatIn, classLabels, C, toler, maxIter,kTup=('lin', 0)):
#full Platt SMO
oS = optStruct(mat(dataMatIn),mat(classLabels).transpose(),C,toler, kTup)
entireSet = T alphaPairsChanged = 0
while (iter & maxIter) and ((alphaPairsChanged & 0) or (entireSet)):
#迭代次数超过指定的最大值或者遍历整个集合都为对任意alpha进行修改时 就退出循环
alphaPairsChanged = 0
if entireSet:
#go over all
遍历所有的值
for i in range(oS.m):
alphaPairsChanged += innerL(i,oS)
#选择第二个alpha并在可能时对其进行优化处理
如果有任意一对alpha值发生变化 返回1
print &fullSet, iter: %d i:%d, pairs changed %d& % (iter,i,alphaPairsChanged)
else:#go over non-bound (railed) alphas
遍历非边界值 不在边界0或C上的值
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)
if entireSet: entireSet = False #toggle entire set loop
elif (alphaPairsChanged == 0): entireSet = True
print &iteration number: %d& % iter
return oS.b,oS.alphas
'''
&&& datMat = mat(dataArr)
&&& datMat[0]*mat(ws)+b
matrix([[-0.]])
如果该值大于0，那么其属于1类;如果该值小于0，那么则属于-1类。
对于数据点0，应该得到的类别标签是-1。可以通过如下命令来确认分类结果的正确性
上述代码只有一个可选的输人参数，该输人参数是高斯径向基函数中的一个用户定义变量。
整个代码主要是由以前定义的函数集合构成的。首先，程序从文件中读人数据集，然后在该数据
集上运行Platt SMO算法，其中核函数的类型为·rbf'o
优化过程结束后，在后面的矩阵数学运算中建立了数据的矩阵副本，并且找出那些非零的
alpha值，从而得到所需要的支持向量;同时，也就得到了这些支持向量和alpha的类别标签值。
这些值仅仅是需要分类的值。
整个代码中最重要的是for循环开始的那两行，它们给出了如何利用核函数进行分类。首先
利用结构初始化方法中使用过的kernelTrans()函数，得到转换后的数据。然后，再用其与前
面的alpha及类别标签值求积。其中需要特别注意的另一件事是，在这几行代码中，是如何做到只
需要支持向量数据就可以进行分类的。除此之外，其他数据都可以直接舍弃。
与第一个for循环相比，第二个for循环仅仅只有数据集不同，后者采用的是测试数据集。
读者可以比较不同的设置在测试集和训练集上表现出的性能。
对数据进行分类处理
&&& ws = svmMLiA.calcWs(alpha,dataArr,labelArr)
array([[ 0.],
&&& datMat = mat(dataArr)
&&& datMat[0]*mat(ws)+b
matrix([[-0.]])
如果该值大于0，那么其属于1；如果该值小鱼0，那么则属于-1
对于数据0 应该得到的类别标签是-1
'''
def calcWs(alphas,dataArr,classLabels):
X = mat(dataArr); labelMat = mat(classLabels).transpose()
m,n = shape(X)
w = zeros((n,1))
for i in range(m):
w += multiply(alphas[i]*labelMat[i],X[i,:].T)
'''
这时观察一下函数testRbf()
的输出结果就会发现，此时的测试错误率也在下降。该数据集在这个设置的某处存在着最优值。
如果降低。，那么训练错误率就会降低，但是测试错误率却会上升。
支持向量的数目存在一个最优值。SVM的优点在于它能对数据进行高效分类。如果支持向量
太少，就可能会得到一个很差的决策边界(下个例子会说明这一点);如果支持向量太多，也就
相当于每次都利用整个数据集进行分类，这种分类方法称为k近邻。
利用核函数进行分类的径向基测试函数
'''
def testRbf(k1=1.3):
#k1 可选参数
高斯径向基函数中的一个用户定义变量
dataArr,labelArr = loadDataSet('testSetRBF.txt')
b,alphas = smoP(dataArr, labelArr, 200, 0., ('rbf', k1)) #C=200 important
datMat=mat(dataArr); labelMat = mat(labelArr).transpose()
svInd=nonzero(alphas.A&0)[0]
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):
kernelEval = kernelTrans(sVs,datMat[i,:],('rbf', k1))
#调用kernelTrans函数进行结构初始化函数，得到转换后的数据
predict=kernelEval.T * multiply(labelSV,alphas[svInd]) + b
#用得到的数据与前面的alpha及类别标签值求积
if sign(predict)!=sign(labelArr[i]): errorCount += 1
print &the training error rate is: %f& % (float(errorCount)/m)
dataArr,labelArr = loadDataSet('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)
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
'''
基于SVM的数字识别：
(1)收集数据:提供的文本文件。
(2)准备数据:基于二值图像构造向量。
(3)分析数据:对图像向量进行目测。
(4)训练算法:采用两种不同的核函数，并对径向基核函数采用不同的设置来运行SMO算法。
(5)测试算法:编写一个函数来测试不同的核函数并计算错误率。
(6)使用算法:一个图像识别的完整应用还需要一些图像处理的知识，这里并不打算深入介绍。
基于SVM的手写数字识别
'''
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
'''
下一个函数七estDigi七。()并不是全新的函数，它和testRbf()的代码几乎一样，唯一的大
区别就是它调用了load工mages()函数来获得类别标签和数据。另一个细小的不同是现在这里的
函数元组kTup是输人参数，而在testRbf()中默认的就是使用rb唯函数。如果对于函数
七estDigi七。()不增加任何输人参数的话，那么kTup的默认值就是('rbf ,10)0
你可能注意到了一个有趣的现象，即最小的训练错误率并不对应于最小的支持向量数目。另
一个值得注意的就是，线性核函数的效果并不是特别的糟糕。可以以牺牲线性核函数的错误率来
换取分类速度的提高。尽管这一点在实际中是可以接受的，但是还得取决于具体的应用。
'''
def testDigits(kTup=('rbf', 10)):
dataArr,labelArr = loadImages('trainingDigits')
b,alphas = smoP(dataArr, labelArr, 200, 0., 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('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)
'''#######********************************
Non-Kernel VErsions below
'''#######******