【学习机器学习】实验——线性模型实现

前言

上周五机器学习的实验课内容,由于图像处理占用了大量时间(即便现在还在研究),所以腾不出时间写博客了。而且机器学习部分的实验内容多数要求调用matlab自带函数,因此也没有什么可讲的编程部分。下面的内容,我直接拷贝了我的实验报告。

一、实验目的

  1. 实现对数几率回归或多元线性回归算法

  2. 实现线性判别LDA算法

二、实验内容

1、导入数据

本次实验我选择调用matlab自带的fisheriris数据集,使用load fisheriris;将其加入工作区,可以看到matlab将其分为了属性meas和标签species两部分:

在这里插入图片描述

并且可以看到标签是以字符串形式保存的,这不利于我们后续的操作,因此我写了一个函数将cell类型的分类标签转换为1、2、3…编号:

function [types]=tag2num(data1)
%将种类名标签,转换为数字,方便与属性合并后传入k折函数
types=[];
data1=categorical(data1);   % 先将cell转换为categorical数组
species=unique(data1);    % 取出类别

for i=1:length(data1)
    for j=1:length(species)
    if data1(i)==species(j)
        types(end+1)=j; % 按照类别给样本标号
    end
    end
end
types=types';   % 转置

再将属性与转换后的标签合并起来,方便我们下一步划分训练集与测试集。

types=tag2num(species);
iris=[meas,types];
2、划分训练集与测试集

这里我并没有选择简单的randperm打乱数据,而是调用了之前写的k折交叉验证函数。因为我们的线性模型,尤其是线性判别分析主要用于二分类问题,因此我将iris数据集中标签值为3的数据删掉了,之后的操作中仅考虑二分类问题。

iris(find(iris(:,5)==3),:)=[];  % 因为我们只处理二分类问题,只保留标记值为12的样本
[Train_kcross,Test_kcross]=kcrossvalidation(iris);  %k折分组

p次k折交叉验证一共会给出100组测试集和训练集,我们只需要选取其中的一组进行后续的调用即可,具体的选取根据函数还有所不同,见后续操作。

3、多元线性回归
3.1 多元线性回归系数矩阵

matlab自带有多元线性回归的实现函数regress。考虑到我们需要可视化最后的结果,而iris本身有4个属性,这意味着我们需要做一个超过三维的图像,显然是不可能的。所以,我们在选取数据时只选用两个属性值作为输入数据。

Train_x=Train_kcross(:,3:4,1,1);   % 将第一组训练样本拆开,注意只选取了属性34
Train_x=[Train_x,ones(length(Train_x),1)];
Train_y=Train_kcross(:,5,1,1);  % x是属性组,y是标记值
Test_x=Test_kcross(:,3:4,1,1);     % 将第一组测试样本拆开
Test_y=Test_kcross(:,5,1,1);    % x是属性组,y是标记值

注意到在输入数据Train_x后我还加了一列1,这取决于你是否希望regress函数输出常数项b,如果没有这一列1则regress默认不输出常数项。关于regress函数,其具体的定义如下:

[B,BINT,R,RINT,STATS] = regress(Y,X)

其返回值的含义是:

B: 回归系数向量。当输入数据X中有一列1时,B(1)将表示常数项(截距);否则,B表示对应属性的系数

BINT:系数估计值的置信度为95%的置信区间

R:残差

RINT:各残差的置信区间

STATS:用于检验回归模型的统计量。有4个数值:判定系数R^2,F统计量观测值,检验的p的值,误差方差的估计。其中检验值p应当小于0.05,否则回归无效。

因此我们分别将Train_y和Train_x作为参数调用regress函数即可。实验中,我调用函数之后得到系数矩阵b=[0.2664,0.1809,0.5959],下面我们要将其可视化。

3.2 多元线性回归可视化

我用三维图像将上面得到的回归模型进行了可视化,可视化代码如下:

% 我们将标签为1和标签为2的分别取出 %
Train_x_class1=Train_kcross(find(Train_kcross(:,5)==1),3:4,1,1);
Train_x_class2=Train_kcross(find(Train_kcross(:,5)==2),3:4,1,1);

% 三维作图 %
figure;
hold on;
z_class1=zeros(length(Train_x_class1),1);
z_class1(:)=1;
z_class2=zeros(length(Train_x_class2),1);
z_class2(:)=2;
scatter3(Train_x_class1(:,1),Train_x_class1(:,2),z_class1,'filled');
scatter3(Train_x_class2(:,1),Train_x_class2(:,2),z_class2,'filled');
x1fit=min(Train_x(:,1)):0.1:max(Train_x(:,1));
x2fit=min(Train_x(:,2)):0.1:max(Train_x(:,2));
[X1FIT,X2FIT]=meshgrid(x1fit,x2fit);
YFIT=b(1)+b(2)*X1FIT+b(3)*X2FIT;
mesh(X1FIT,X2FIT,YFIT);
view(10,10);
xlabel('x1'); %设置x轴的名称  
ylabel('x2'); %设置y轴的名称  
zlabel('y');  %设置z轴的名称
hold off;

图像如下:

多元线性回归可视化
实际上划分的效果比较一般,这张图已经是通过调整角度使它“看起来”比较成功了。

3.3 测试集检验

从上面的图像看来,我们得到的回归模型分类结果理应有一定的误差,我们选用一些测试集来计算一下精度。

regress_predict=b(2)*Test_x(:,1)+b(3)*Test_x(:,2)+b(1);	% 计算预测结果 

这里我不清楚如何选取判定用的阈值,如果把1作为阈值,那么精度往往是100%。我问过老师,老师说iris数据集属于“太好分”了,所以对模型要求不高,大家可以试试用别的数据集。

% 计算精度和错误率
right=0;
wrong=0;
for i=1:length(regress_predict)
    if(regress_predict(i)<=1)
        regress_predict(i)=1;
    end
    if(regress_predict(i)>1 && regress_predict(i)<3)
        regress_predict(i)=2;
    end
end
for i=1:length(regress_predict)
    if regress_predict(i)==Test_y(i)
        right=right+1;
    else
        wrong=wrong+1;
    end
end

error_rate=wrong/length(regress_predict);
accuracy=right/length(regress_predict);
4、线性判别分析
4.1 编写线性判别分析

我并不确定matlab自带函数classify()是不是用的我们所讲的线性判别分析的方法,因此选择自己写了一个。LDA算法的根本就在于最大化目标J:

J = ∣ ∣ ω T μ 0 − ω T μ 1 ∣ ∣ 2 ω T ∑ 0 ω + ω T ∑ 1 ω = ω T ( μ 0 − μ 1 ) ( μ 0 − μ 1 ) T ω ω T ( ∑ 0 + ∑ 1 ) ω J=\frac{||\omega^{T}\mu_{0} -\omega^{T}\mu_{1}||^{2}}{\omega^{T}\sum_{0}\omega + \omega^{T}\sum_{1}\omega}=\frac{\omega^{T}(\mu_{0}-\mu_{1})(\mu_{0}-\mu_{1})^{T}\omega}{\omega^{T}(\sum_{0}+\sum_{1})\omega} J=ωT0ω+ωT1ωωTμ0ωTμ12=ωT(0+1)ωωT(μ0μ1)(μ0μ1)Tω

我们将 ∑ 0 \sum_{0} 0 ∑ 1 \sum_{1} 1分别记为s1和s2,那么就有类内散度矩阵 S ω = s 1 + s 2 S_{\omega}=s1+s2 Sω=s1+s2,经过一系列推到之后我们知道结论,我们所求的 ω = S ω − 1 ( μ 0 − μ 1 ) \omega=S_{\omega}^{-1}(\mu_{0}-\mu_{1}) ω=Sω1(μ0μ1),我们的函数编写也是相同的思路,具体代码如下:

function [output_class,w] = LDA(input_sample, class1_samples, class2_samples)
% 形参解释 %
% input_sample 测试样例,维度为N
% class1_sample 类别一的训练样例,样本个数NC1,维度为N
% class2_sample 类别二的训练样例,样本个数NC2,维度为N
% output_class 返回测试样例的预测结果,应为1/2 %
% w 权值/系数矩阵 %
[m1, n1] = size(class1_samples);
[m2, n2] = size(class2_samples);
mu1 = sum(class1_samples) / size(class1_samples, 1);    % 均值μ1
mu2 = sum(class2_samples) / size(class2_samples, 1);    % 均值μ2
s1 = 0;
s2 = 0;

for i = 1:m1
    s1 = s1 + (class1_samples(i,:) - mu1)' * (class1_samples(i,:) - mu1);
end
for i = 1:m2
    s2 = s2 + (class2_samples(i,:) - mu2)' * (class2_samples(i,:) - mu2);
end
sw = s1 + s2;   % 类内散度矩阵
w = inv(sw)' * (mu1 - mu2)';    % 结果w(系数矩阵)
separationPoint = (mu1 + mu2) * w * 0.5;    % 分隔点为1/2 * w * (μ1+μ2)
output_class = (input_sample * w < separationPoint) + 1;
% bool的结果为01,因此这里得到的是分类值12

调用这个函数[output_class,w] = LDA(Test_x, Train_x_class1, Train_x_class2);就可以得到投影方向矩阵w和测试样例的预测结果output_class。

4.2 线性判别分析可视化

线性判别我们同样是选取两个属性来做图,而且,通过不同属性的作图我们实际上可以区分出一个属性是“好”还是“坏”。代码如下:

% 下面选属性作图,34是比较好的属性 %
figure;
hold on;
plot(Train_x_class1(:,3),Train_x_class1(:,4),'r+','markerfacecolor',[1,0,0]);
plot(Train_x_class2(:,3),Train_x_class2(:,4),'b*','markerfacecolor',[0,0,1]);
grid on;
k2=w(3)/w(4);
x2=1:6;
y2=k2*x2;
plot(x2,y2);
hold off;

% 12是比较差的属性 %
figure;
hold on;
plot(Train_x_class1(:,1),Train_x_class1(:,2),'r+','markerfacecolor',[1,0,0]);
plot(Train_x_class2(:,1),Train_x_class2(:,2),'b*','markerfacecolor',[0,0,1]);
grid on;
k3=w(1)/w(2);
x3=1:8;
y3=k3*x3;
plot(x3,y3);
hold off;

我们可以看一下效果:

在这里插入图片描述

左图是使用1、2属性绘制的图像,可以看到点到直线的投影混杂在一起,并没有达到区分的目的,这说明1、2两个属性并不好。我举一个可能不太恰当的例子,1、2属性就像是性别和身高,而输出结果是判断一个人是好人坏人,由于这个属性本身对于结果的影响很小,因此用其作为判别依据并不好。

而右图是使用3、4属性绘制的图像,显然是很明显的区分了两个类,这样的属性就是好属性。

4.3 测试集检验

我写的函数在调用时就输出了测试样例的预测结果,我们按照与多元线性回归相同的方法统计精度和错误率。

在这里插入图片描述

精度达到了100%,并且我多次调用结果均为如此,不知是确实有如此高的精度还是在统计过程中有些问题。

三、实验总结

本次实验我们完成了多元线性回归和线性判别分析两种线性模型的算法。线性模型常用于处理分类问题,如果我们的统计方法没有问题,那么可以看到这两种线性模型在进行分类时都有非常高的精度,而且更重要的是这两种算法的实现并不复杂,因此线性模型现在依旧是机器学习中一个重要的部分。

PS:对数几率回归

由于实验限两天内上交,所以当时没有写对数几率回归。后来抽时间实现了一下,参考的是这篇博客

% logistic regression %
load fisheriris;    % 用鸢尾花数据集,得到meas和species
types=tag2num(species); % 把species中的标签名转换为123
iris=[meas,types];     % 把属性与标签合并
iris(find(iris(:,5)==3),:)=[];  % 因为我们只处理二分类问题,只保留标记值为12的样本
[Train_kcross,Test_kcross]=kcrossvalidation(iris);  %k折分组

% 选出一组训练集和测试集 %
Train_x=Train_kcross(:,3:4,1,1);   % 将第一组训练样本拆开
Train_y=Train_kcross(:,5,1,1);  % x是属性组,y是标记值
Test_x=Test_kcross(:,3:4,1,1);     % 将第一组测试样本拆开
Test_y=Test_kcross(:,5,1,1);    % x是属性组,y是标记值

% 我们将标签为1和标签为2的分别取出 %
Train_x_class1=Train_kcross(find(Train_kcross(:,5)==1),3:4,1,1);
Train_x_class2=Train_kcross(find(Train_kcross(:,5)==2),3:4,1,1);

figure;
hold on;
plot(Train_x_class1(:,1),Train_x_class1(:,2),'k+','LineWidth',2,'MarkerSize',7);
plot(Train_x_class2(:,1),Train_x_class2(:,2),'ko','MarkerFaceColor','y','MarkerSize',7);
xlabel('x1');
ylabel('x2');
hold off;

[m,n]=size(Train_x);
Train_x=[ones(m,1) Train_x];
initial_theta = zeros(n + 1, 1);
options = optimset('GradObj', 'on', 'MaxIter', 400);
[theta, cost] = fminunc(@(t)(costFunction(t, Train_x, Train_y)), initial_theta, options);

figure;
hold on;
if size(Train_x, 2) <= 3
    % Only need 2 points to define a line, so choose two endpoints
    plot_x = [min(Train_x(:,2))-2,  max(Train_x(:,2))+2];
    % Calculate the decision boundary line
    plot_y = (-1./theta(3)).*(theta(2).*plot_x + theta(1));
    plot(plot_x, plot_y)
    % Legend, specific for the exercise
    legend('Admitted', 'Not admitted', 'Decision Boundary')
    axis([30, 100, 30, 100])
else
    % Here is the grid range
    u = linspace(-1, 1.5, 50);
    v = linspace(-1, 1.5, 50);
    z = zeros(length(u), length(v));
    % Evaluate z = theta*x over the grid
    for i = 1:length(u)
        for j = 1:length(v)
            z(i,j) = mapFeature(u(i), v(j))*theta;
        end
    end
    z = z'; 
    contour(u, v, z, [0, 0], 'LineWidth', 2)
end
hold off;

其中的costFunction定义如下

%costFunction.m
function [J, grad] = costFunction(theta, X, y)
m = length(y); % number of training examples
J = 0;
grad = zeros(size(theta));
h=1.0./(1.0+exp(-1*X*theta));
m=size(y,1);
J=((-1*y)'*log(h)-(1-y)'*log(1-h))/m;
for i=1:size(theta,1)
    grad(i)=((h-y)'*X(:,i))/m;
end
  • 1
    点赞
  • 25
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值