机器学习实战

来源：互联网发布：mt4顾比均线指标源码编辑：程序博客网时间：2024/06/05 20:53

机器学习实战ByMatlab（1）：KNN算法

KNN 算法其实简单的说就是“物以类聚”，也就是将新的没有被分类的点分类为周围的点中大多数属于的类。它采用测量不同特征值之间的距离方法进行分类，思想很简单：如果一个样本的特征空间中最为临近（欧式距离进行判断）的K个点大都属于某一个类，那么该样本就属于这个类。这就是物以类聚的思想。

当然，实际中，不同的K取值会影响到分类效果，并且在K个临近点的选择中，都不加意外的认为这K个点都是已经分类好的了，否则该算法也就失去了物以类聚的意义了。

KNN算法的不足点：

1、当样本不平衡时，比如一个类的样本容量很大，其他类的样本容量很小，输入一个样本的时候，K个临近值中大多数都是大样本容量的那个类，这时可能就会导致分类错误。改进方法是对K临近点进行加权，也就是距离近的点的权值大，距离远的点权值小。

2、计算量较大，每个待分类的样本都要计算它到全部点的距离，根据距离排序才能求得K个临近点，改进方法是：先对已知样本点进行剪辑，事先去除对分类作用不大的样本。

适用性：

适用于样本容量比较大的类域的自动分类，而样本容量较小的类域则容易误分

算法描述：

1、计算已知类别数据集合汇总的点与当前点的距离
2、按照距离递增次序排序
3、选取与当前点距离最近的K个点
4、确定距离最近的前K个点所在类别的出现频率
5、返回距离最近的前K个点中频率最高的类别作为当前点的预测分类

Matlab 实现

这里以一个完整实例呈现，代码如下：

function relustLabel = KNN(inx,data,labels,k)

% inx 为输入测试数据，data为样本数据，labels为样本标签

[datarow , datacol] = size(data);

diffMat = repmat(inx,[datarow,1]) - data ;

distanceMat = sqrt(sum(diffMat.^2,2));

[B , IX] = sort(distanceMat,'ascend');

len = min(k,length(B));

relustLabel = mode(labels(IX(1:len)));

end

可以看到，整个KNN算法的Matlab代码也就只有6行，比Python少很多，这其中要得益于Matlab成熟的矩阵计算和很多成熟的函数。

实际调用中，我们利用一个数据集进行测试，该数据集是由1000个样本的3维坐标组成，共分为3个类

首先可视化我们的数据集，看看它的分布：

第一维和第二维：可以清晰地看到数据大致上分为 3 类

第1维和第3维：从这个角度看，3类的分布就有点重叠了，这是因为我们的视角原因

画出3维，看看它的庐山真面目：

由于我们已经编写了KNN代码，接下来我们只需调用就行。了解过机器学习的人都应该知道，很多样本数据在代入算法之前都应该进行归一化，这里我们将数据归一化在[0,1]的区间内，归一化方式如下：

newData =(oldData-minValue)/(maxValue-minValue)

其中，maxValue为oldData的最大值，minValue为 oldData 的最小值。

同时，我们对于1000个数据集，采取10%的数据进行测试，90%的数据进行训练的方式，由于本测试数据之间完全独立，可以随机抽取10%的数据作为测试数据，代码如下：

257

function KNNdatgingTest

clc

clear

close all

data = load('datingTestSet2.txt');

dataMat = data(:,1:3);

labels = data(:,4);

len = size(dataMat,1);

k = 4;

error = 0;

% 测试数据比例

Ratio = 0.1;

numTest = Ratio * len;

% 归一化处理

maxV = max(dataMat);

minV = min(dataMat);

range = maxV-minV;

newdataMat = (dataMat-repmat(minV,[len,1]))./(repmat(range,[len,1]));

% 测试

for i = 1:numTest

classifyresult = KNN(newdataMat(i,:),newdataMat(numTest:len,:),labels(numTest:len,:),k);

fprintf('测试结果为：%d 真实结果为：%d\n',[classifyresult labels(i)])

if(classifyresult~=labels(i))

error = error+1;

end

fprintf('准确率为：%f\n',1-error/(numTest))

end

当我们选择K为4的时候，准确率为：97%

KNN进阶

接下来我们将运用KNN算法实现一个手写识别系统，训练数据集大约2000个样本，每个数字大概有200个样本
测试数据大概有900个样本，由于每个样本都是一个32×32的数字，我们将其转换为1×1024的矩阵，方便我们利用KNN算法
数据如下：

由于数据量比较大，加载数据的时候回花一点时间，具体代码如下：

function handWritingTest

clc

clear

close all

%% 获取目录下的所有txt文件名称

d = dir(['digits/trainingDigits/' '*.txt']); % struct 类型

dircell = struct2cell(d); %cell 类型

trainSetLen = size(dircell,2);

K = 4;

dataSize = 1024;

trainLabels = zeros(trainSetLen,1);

trainSet = [];

simpleTrainSet = zeros(1,dataSize);

simpleTestSet = zeros(1,dataSize);

%% 加载数据

fprintf('loading data...')

for i = 1:trainSetLen

trainName = dircell(1,i);

trainFilename = cell2mat(trainName);

trainLabels(i) = str2num(trainFilename(1));

fid = fopen(['digits/trainingDigits/' trainFilename],'r');

traindata = fscanf(fid,'%s');

for j = 1:dataSize

simpleTrainSet(j) = str2num(traindata(j));

end

trainSet = [trainSet ; simpleTrainSet];

fclose(fid);

end

d = dir(['digits/testDigits/' '*.txt']); % struct 类型

dircell = struct2cell(d); %cell 类型

testSetLen = size(dircell,2);

error = 0;

%% 测试数据

for k = 1:testSetLen

testName = dircell(1,k);

testFilename = cell2mat(testName);

testLabels = str2num(testFilename(1));

fid = fopen(['digits/testDigits/' testFilename],'r');

testdata = fscanf(fid,'%s');

for j = 1:dataSize

simpleTestSet(j) = str2num(testdata(j));

end

classifyResult = KNN(simpleTestSet,trainSet,trainLabels,K);

fprintf('识别数字为：%d 真实数字为：%d\n' , [classifyResult , testLabels])

if(classifyResult~=testLabels)

error = error+1;

end

fclose(fid);

end

fprintf('识别准确率为：%f\n',1-error/testSetLen)

end

不同的K识别准确率稍有不同，当K为4的时候，准确率为 98.31%

机器学习实战ByMatlab（2）：PCA算法

PCA 算法也叫主成分分析（principal components analysis），主要是用于数据降维的。

为什么要进行数据降维？因为实际情况中我们的训练数据会存在特征过多或者是特征累赘的问题，比如：

· 一个关于汽车的样本数据，一个特征是”km/h的最大速度特征“，另一个是”英里每小时“的最大速度特征，很显然这两个特征具有很强的相关性

· 拿到一个样本，特征非常多，样本缺很少，这样的数据用回归去你和将非常困难，很容易导致过度拟合

PCA算法就是用来解决这种问题的，其核心思想就是将 n 维特征映射到 k 维上（k< n），这 k 维是全新的正交特征。我们将这 k 维成为主元，是重新构造出来的 k 维特征，而不是简单地从 n 维特征中取出其余 n-k 维特征。

PCA 的计算过程

假设我们得到 2 维数据如下：

其中行代表样例，列代表特征，这里有10个样例，每个样例有2个特征，我们假设这两个特征是具有较强的相关性，需要我们对其进行降维的。

第一步：分别求 x 和 y 的平均值，然后对所有的样例都减去对应的均值

这里求得 x 的均值为 1.81 ， y 的均值为 1.91，减去均值后得到数据如下：

注意，此时我们一般应该在对特征进行方差归一化，目的是让每个特征的权重都一样，但是由于我们的数据的值都比较接近，所以归一化这步可以忽略不做

第一步的算法步骤如下：

本例中步骤3、4没有做。

第二步：求特征协方差矩阵

公式如下：

第三步：求解协方差矩阵的特征值和特征向量

第四步：将特征值从大到小进行排序，选择其中最大的 k 个，然后将其对应的 k 个特征向量分别作为列向量组成特征矩阵

这里的特征值只有两个，我们选择最大的那个，为：1.28402771 ，其对应的特征向量为：

注意：matlab 的 eig 函数求解协方差矩阵的时候，返回的特征值是一个特征值分布在对角线的对角矩阵，第i 个特征值对应于第 i 列的特征向量

第五步：将样本点投影到选取的特征向量上

假设样本列数为 m ，特征数为 n ，减去均值后的样本矩阵为 DataAdjust(m*n),协方差矩阵为 n*n ,选取 k 个特征向量组成后的矩阵为 EigenVectors(n*k)，则投影后的数据 FinalData 为：

FinalData （m*k） = DataAdjust(m*n) X EigenVectors(n*k)

得到的结果是：

这样，我们就将 n 维特征降成了 k 维，这 k 维就是原始特征在 k 维上的投影。

整个PCA的过程貌似很简单，就是求协方差的特征值和特征向量，然后做数据转换。但为什么协方差的特征向量就是最理想的 k 维向量？这个问题由PCA的理论基础来解释。

PCA 的理论基础

关于为什么协方差的特征向量就是 k 维理想特征，有3个理论，分别是：

1. 最大方差理论

2. 最小错误理论

3. 坐标轴相关度理论

这里简单描述下最大方差理论：

最大方差理论

信号处理中认为信号具有较大的方差，噪声有较小的方差，信噪比就是信号与噪声的方差比，越大越好。因此我们认为，最好的 k 为特征既是将 n 维样本点转换为k 维后，每一维上的样本方差都很大

PCA 处理图解如下：

降维转换后：

上图中的直线就是我们选取的特征向量，上面实例中PCA的过程就是将空间的2维的点投影到直线上。

那么问题来了，两幅图都是PCA的结果，哪一幅图比较好呢？

根据最大方差理论，答案是左边的图，其实也就是样本投影后间隔较大，容易区分。

其实从另一个角度看，左边的图每个点直线上的距离绝对值之和比右边的每个点到直线距离绝对值之和小，是不是有点曲线回归的感觉？其实从这个角度看，这就是最小误差理论：选择投影后误差最小的直线。

再回到上面的左图，也就是我们要求的最佳的 u ，前面说了，最佳的 u 也就是最佳的曲线，它能够使投影后的样本方差最大或者是误差最小。

另外，由于我们前面PCA算法第一步的时候已经执行对样本数据的每一维求均值，并让每个数据减去均值的预处理了，所以每个特征现在的均值都为0，投影到特征向量上后，均值也为0.因此方差为：

最后的等式中中间的那部分其实就是样本方差的协方差矩阵（xi 的均值为 0）

由于 u 是单位向量，得到

上式两边痛乘以 u，得到：

于是我们得到

最佳投影直线就是特征值 λ 最大是对应的特征向量，其次是 λ 第二大对应的特征向量（求解的到的特征向量都是正交的）。其中λ 就是我们的方差，也对应了我们前面的最大方差理论，也就是找到能够使投影后方差最大的直线。

Matlab 实现

function [lowData,reconMat] = PCA(data,K)

[row , col] = size(data);

meanValue = mean(data);

%varData = var(data,1,1);

normData = data - repmat(meanValue,[row,1]);

covMat = cov(normData(:,1),normData(:,2));%求取协方差矩阵

[eigVect,eigVal] = eig(covMat);%求取特征值和特征向量

[sortMat, sortIX] = sort(eigVal,'descend');

[B,IX] = sort(sortMat(1,:),'descend');

len = min(K,length(IX));

eigVect(:,IX(1:1:len));

lowData = normData * eigVect(:,IX(1:1:len));

reconMat = (lowData * eigVect(:,IX(1:1:len))') + repmat(meanValue,[row,1]); % 将降维后的数据转换到新空间

end

调用方式

function testPCA

clc

clear

close all

filename = 'testSet.txt';

K = 1;

data = load(filename);

[lowData,reconMat] = PCA(data,K);

figure

scatter(data(:,1),data(:,2),5,'r')

hold on

scatter(reconMat(:,1),reconMat(:,2),5)

hold off

end

效果图

机器学习实战ByMatlab（5）：LogisticRegression

什么叫做回归呢？举个例子，我们现在有一些数据点，然后我们打算用一条直线来对这些点进行拟合（该曲线称为最佳拟合曲线），这个拟合过程就被称为回归。

利用Logistic回归进行分类的主要思想是：

根据现有数据对分类边界线建立回归公式，以此进行分类。

这里的”回归“一词源于最佳拟合，表示要找到最佳拟合参数集。训练分类器时的嘴阀就是寻找最佳拟合曲线，使用的是最优化算法。

基于Logistic回归和Sigmoid函数的分类

优点：计算代价不高，易于理解和实现
缺点：容易欠拟合，分类精度可能不高

使用数据类型：数值型和标称型数据

Sigmoid函数：

波形如下：

当z为0时，值为0.5，当z增大时，g(z)逼近1，当z减小时，g(z)逼近0

Logistic回归分类器：

对每一个特征都乘以一个回归系数，然后把所有结果都相加，再讲这个总和代入Sigmoid函数中，从而得到一个范围在0-1之间的数值。任何大于0.5的数据被分为1，小于0.5的数据被分为0.因此Logistic回归也被看成是一种概率分布。

分类器的函数形式确定之后，现在的问题就是，如何确定回归系数？

基于最优化方法的最佳回归系数确定

Sigmoid函数的输入记为z，由下面公式得出：

如果采用向量的写法，则上述公式可以写成：

其中向量X就是分类器的输入数据，向量W也就是我们要找到的最佳参数，从而使分类器尽可能更加地精确。接下来将介绍几种需找最佳参数的方法。

梯度上升法

梯度上升法的基本思想：

要找到某函数的最大值，最好的方法是沿着该函数的梯度方向寻找

这里提一下梯度下降法，这个我们应该会更加熟悉，因为我们在很多代价函数J的优化的时候经常用到它，其基本思想是：

要找到某函数的最小值，最好的方法是沿着该函数的梯度方向的反方向寻找

函数的梯度表示方法如下：

移动方向确定了，移动的大小我们称之为步长，用α表示，用向量来表示的话，梯度下降算法的迭代公式如下：

该公式已知被迭代执行，直到某个停止条件位置，比如迭代次数达到某个指定值或者算法的误差小到某个允许的误差范围内。

注：梯度下降算法中的迭代公式如下：

Matlab实现

function weight = gradAscent

clc

close all

clear

data = load('testSet.txt');

[row , col] = size(data);

dataMat = data(:,1:col-1);

dataMat = [ones(row,1) dataMat] ;

labelMat = data(:,col);

alpha = 0.001;

maxCycle = 500;

weight = ones(col,1);

for i = 1:maxCycle

h = sigmoid((dataMat * weight)');

error = (labelMat - h');

weight = weight + alpha * dataMat' * error;

end

figure

scatter(dataMat(find(labelMat(:) == 0),2),dataMat(find(labelMat(:) == 0),3),3);

hold on

scatter(dataMat(find(labelMat(:) == 1),2),dataMat(find(labelMat(:) == 1),3),5);

hold on

x = -3:0.1:3;

y = (-weight(1)-weight(2)*x)/weight(3);

plot(x,y)

hold off

end

function returnVals = sigmoid(inX)

% 注意这里的sigmoid函数要用点除

returnVals = 1.0./(1.0+exp(-inX));

end

效图如下：

由上图可以看到，回归效果还是挺不错的，只有2-4个点分类错误。

其实这是的梯度上升算法是批量梯度上升算法，每一次更新参数的时候都要讲所有的数据集都代入训练，效果并不好，下面我们将介绍改进版本：随机梯度上升算法

随机梯度上升

梯度上升算法在每次更新回归系数时都要遍历整个数据集，该方法在处理100个左右的数据集时尚可，但如果有数十亿样本和成千上万的特征，那么该方法的复杂度就太高了。一种改进方法是一次仅用一个样本点来更新回归系数，该方法就称为随机梯度上升法。由于可以在新样本到来之前对分类器进行增量式更新，因此随机梯度算法是一个在线学习算法。与”在线学习“相对应，一次处理所有数据被称作是”批处理“

随机梯度上升算法可以写成如下的伪代码：

所有回归系数初始化为1
对数据集中的每个样本
计算该样本的梯度
使用alpha x gradient 更新回归系数值
返回回归系数值

Matlab 代码实现

function stocGradAscent

% Description : LogisticRegression using stocGradAsscent

% Author : Liulongpo

% Time：2015-4-18 10:57:25

clc

clear

close all

data = load('testSet.txt');

[row , col] = size(data);

dataMat = [ones(row,1) data(:,1:col-1)];

alpha = 0.01;

labelMat = data(:,col);

weight = ones(col,1);

for i = 1:row

h = sigmoid(dataMat(i,:)*weight);

error = labelMat(i) - h;

dataMat(i,:)

weight

weight = weight + alpha * error * dataMat(i,:)'

end

figure

scatter(dataMat(find(labelMat(:)==0),2),dataMat(find(labelMat(:)==0),3),5);

hold on

scatter(dataMat(find(labelMat(:) == 1),2),dataMat(find(labelMat(:) == 1),3),5);

hold on

x = -3:0.1:3;

y = -(weight(1)+weight(2)*x)/weight(3);

plot(x,y)

hold off

end

function returnVals = sigmoid(inX)

% 注意这里的sigmoid函数要用点除

returnVals = 1.0./(1.0+exp(-inX));

end

效果如下：

由上图可以看出，随机梯度上升算法分类效果并没有上面的的梯度上升算法分类效果好。

但是直接比较梯度上升算法和随机梯度上升算法是不公平的，前者是在整个数据集上迭代500次得到的结果，后者只是迭代了100次。一个判断算法优劣的可靠方法是看它是否收敛，也就是说求解的参数是否达到了稳定值，是否还会不断变化。

我们让随机梯度上升算法在整个数据集上运行200次，迭代过程中3个参数的变化如下图：

由上图可以看到，weight1 最先达到稳定，而weight0和weight2则还需要更多的迭代次数来达到稳定。

此时的分类器跟之前的梯度上升算法的分类效果差不多，如下：

但同时我们也可以看到，三个参数都有不同程度的波动。产生这种现象的原因是存在一些不能被正确分类的样本点（数据集并非线性可分），在每次迭代的时候都会引起参数的剧烈变化。我们期望算法能避免来回波动，从而收敛到某个值。另外，算法收敛速度也要加快。

改进的随机梯度上升算法

改进的随机梯度上升算法的主要两个改进点如下：

1,每一步调整alpha的值，也就是alpha的值是不严格下降的
2.随机采取样本来更新回归参数

matlab代码如下：

function ImproveStocGradAscent

% Description : LogisticRegression using stocGradAsscent

% Author : Liulongpo

% Time：2015-4-18 10:57:25

clc

clear

close all

data = load('testSet.txt');

[row , col] = size(data);

dataMat = [ones(row,1) data(:,1:col-1)];

%alpha = 0.01;

numIter = 20;

labelMat = data(:,col);

weightVal = zeros(3,numIter*row);

weight = ones(col,1);

j = 0;

for k = 1:numIter

randIndex = randperm(row);

for i = 1:row

% 改进点 1

alpha = 4/(1.0+i+k)+0.01;

j = j+1;

% 改进点 2

h = sigmoid(dataMat(randIndex(i),:)*weight);

% 改进点 2

error = labelMat(randIndex(i)) - h;

% 改进点 2

weight = weight + alpha * error * dataMat(randIndex(i),:)';

weightVal(1,j) = weight(1);

weightVal(2,j) = weight(2);

weightVal(3,j) = weight(3);

end

figure

i = 1:numIter*row;

subplot(3,1,1)

plot(i,weightVal(1,:)),title('weight0')%,axis([0 numIter*row 0.8 7])

j = 1:numIter*row;

subplot(3,1,2)

plot(j,weightVal(2,:)),title('weight1')%,axis([0 numIter*row 0.3 1.2])

k = 1:numIter*row;

subplot(3,1,3)

plot(k,weightVal(3,:)),title('weight2')%,axis([0 numIter*row -1.2 -0.1])

figure

scatter(dataMat(find(labelMat(:)==0),2),dataMat(find(labelMat(:)==0),3),5);

hold on

scatter(dataMat(find(labelMat(:) == 1),2),dataMat(find(labelMat(:) == 1),3),5);

hold on

x = -3:0.1:3;

y = -(weight(1)+weight(2)*x)/weight(3);

plot(x,y,'r')

hold off

end

function returnVals = sigmoid(inX)

% 注意这里的sigmoid函数要用点除

returnVals = 1.0./(1.0+exp(-inX));

end

改进点 1 中的alpha会随着迭代次数的增加不断减小，但由于代码中常数0.01的存在，alpha不会减少到0。这样做是为了保证在多次迭代之后新数据对于参数的更新还有一定的影响。

另一点值得注意的就是，alpha每次减少 1/(k+i) ，k 是迭代次数，i是样本的下标。所以 alpha 不是严格下降的。避免参数的严格下降也常见于模拟退火算法等其他优化算法中。

第二个改进的地方如代码注释中标记的，这里通过随机采取样本来更新回归参数，这样能够减少参数的周期性的波动。

由于alpha的动态变化，我们可以在开始的时候设置比较大的值，代码中设置0.01，alpha也就是每一次迭代的步长，步长越大，越能够加快参数的收敛速度。然后ahpha会不严格下降，这样就避免了过拟合现象的发生。至于什么是过拟合已经alpha的选取问题将在下面描述。

迭代20次后效果如下：

由上图可知，步长alpha动态变化之后，参数的收敛速度加快了很多，这里只是对所有样本数据集迭代20次，weight0 和 weight2很早就收敛。证明了该算法的优异性。

学习率alpha的选取

首先我们看一下梯度上升算法的核心代码，如下：

h =sigmoid(dataMat(i,:) * weight);
error = labelMat(i) – h;
weight = weight + alpha * error * dataMat(i,:)’;

第一行做的就是估计分类，第二行计算当前估计与正确分类之间的差error，第三行根据这个error来更新参数weight。

我们在迭代的时候，要做的目标就是最小化 error，我们令 J 代表 error，令向量 θ 代表weight，则很显然，J是θ的函数。这里盗用Standfor 机器学习教程的图，如下：

上图中的每个箭头就是每一次迭代的更新步长，第一幅图我们看到，在最小化 J(θ) 的时候迭代了很多次，这说明什么？说明我们要走很多步才能到达全局最优点，原因就是我们每一步走的距离太短，走得太慢，也就是我们的alpha设置得太小。但是当我们处于最优点附近的时候，这样有利我们向最优点靠近。

下图中的每个箭头也代表走一步，我们可以看到，迭代的时候，每一步都没有到达最优点，而是在最优点的附近波动。为什么呢？因为步长太大了嘛，明明就在眼前了，半步或者四分之三步就走到了，你却只能一跨而过，重新再来。但是学习率大的话，在刚开始迭代的时候有利于我们参数的快速收敛，也有利于我们避开局部最小值。

综合以上两种情况，我们就应该在开始的时候选取较大的学习率，然后不断不严格减小学习率，这样才是最优的选择。

那么，我们开始的学习率应该怎么选取？Andrew Ng 在课程中建议先试试0.01，太大就0.003，太小就0.03….

0 0