吴恩达机器学习笔记（2）Logistic回归

这一部分是对吴恩达机器学习logistic回归部分内容的总结，主要分为以下几个部分

1.数据绘制

2.二元Logistic回归

3.相关结论的可视化

4.多项式Logistic回归

1.数据绘制

由于一些基本操作在第一讲都有详细的描述，这一节主要以代码和注释为主

data = load('ex2data1.txt');%读取数据
X = data(:, [1, 2]); y = data(:, 3);%前两个变量为输入，后一个变量为输出
%注意：logistic回归输出只有0和1

这种点状图的画法比较复杂，所以专门定义一个函数

function plotData(X, y)
figure; %建立一个新画布
hold on;%锁定画布，接下来的绘图不会产生新图片
pos = find(y == 1); neg = find(y == 0);
%find函数是一个很重要且很常见的函数，这里作特别说明
%find()函数的功能是返回括号数组或矩阵中满足条件的值所在的位置
%例如X=[1,3,5,7,9]
%find(X==3)
%则会输出2
%通过这个函数我们就能有效地找到y=0的样本位置和y=1的样本位置
%并将得到的位置储存起来方便再次调用
plot(X(pos, 1), X(pos, 2), 'k+','LineWidth', 2, 'MarkerSize', 7);
%X(pos,1)的意思是pos行，第一列，pos就是刚才找到的y=1的值，这种调用很方便
%X(pos,2)同理
%‘k+’是形状，这里表示十字
%‘linewidth是线宽，即十字线宽’ 是2
%‘markersize’是点的大小 是7
plot(X(neg, 1), X(neg, 2), 'ko', 'MarkerFaceColor', 'y','MarkerSize', 7)
%‘Makerfacecolor’即圆形里面的颜色是 y=yellow
hold off;%解锁画布
end

find函数更详细的解释见
https://www.cnblogs.com/anzhiwu815/p/5907033.html
这样，我们就定义了一个绘图函数
则主程序为

plotData(X, y);
hold on;%锁定画布，进行必要的标注
xlabel('Exam 1 score') %x轴
ylabel('Exam 2 score') %y轴
legend('Admitted', 'Not admitted') %图例
hold off;%解锁

这样我们就完成了基本的图形绘制，效果如下

2.二元Logistic回归

有了上述预处理，接下来就可以进行Logistic回归了，所谓Logistic回归，就是给定一些参数，通过算法，来判断基于这种参数的输出是0还是1
算法步骤依然是

1.作出假设函数

具体的公式见笔记，这里仅讨论其代码实现
首先定义西格蒙函数，即假设函数

function g = sigmoid(z)
%其中输入z就是theta*X
g = 1./(1 + exp(-z));
%这里说明，但凡除数和被除数里面有数组，那一定用点除，因为矩阵没有除法
end

2.计算代价函数与梯度

有了假设函数，就可以计算代价函数了，至于为什么梯度也要在这里计算，是因为这次要用一种fminunc的特殊用法，这种用法很有意义
代价函数公式见笔记，这里仅讨论其代码实现
定义这样一个有两个返回值的函数

function [J, grad] = costFunction(theta, X, y)
m = length(y); %样本容量
J = -1 * sum( y .* log( sigmoid(X*theta) ) ...
    + (1 - y ) .* log( (1 - sigmoid(X*theta)) ) ) / m ;
%计算代价函数
grad = ( X' * (sigmoid(X*theta) - y ) )/ m ;
%这里所谓的梯度仅指代价函数J(θ)对θ的偏导数，而非θ的迭代
end

有了这样的函数，主函数就是

[m, n] = size(X);%样本容量的初始化
X = [ones(m, 1) X];%添加x_0便于后面的计算
initial_theta = zeros(n + 1, 1);%theta的初始化
[cost, grad] = costFunction(initial_theta, X, y);
%调用函数计算代价函数和梯度

3.梯度下降最小化代价函数

使得代价函数最小的theta就是最终的估计模型
这里详细介绍一种fminunc的高级用法

options = optimset('GradObj', 'on', 'MaxIter', 400);
%这种用法首先要设置options，其中'GradObj' 'on'代表用户自己选择梯度
%'MaxIter'，'400'就是最大迭代次数为400
[theta, cost] = ...
	fminunc(@(t)(costFunction(t, X, y)), initial_theta, options);
%调用fminunc,其中(@t)就是说t是自变量，只改变t使得第一个参数即costFunction最小
%initial_theta是它的初始值
%最后一个参数表示加载之前的设置

这种用法不需要我们自己设置theta的迭代，我们只要自行这是梯度就能算出theta的最优解
当然，算出theta之后我们要做的就是输出、可视化、测试、计算误差
输出当然就是判断西格蒙函数是否大于0.5，定义函数如下：

function p = predict(theta, X)
%PREDICT Predict whether the label is 0 or 1 using learned logistic 
%regression parameters theta
%   p = PREDICT(theta, X) computes the predictions for X using a 
%   threshold at 0.5 (i.e., if sigmoid(theta'*x) >= 0.5, predict 1)

m = size(X, 1); % Number of training examples

% You need to return the following variables correctly
p = zeros(m, 1);%先全初始化为0
k = find(sigmoid( X * theta) >= 0.5 );%如果西格蒙函数>=0.5
p(k)= 1;%则对应的输出为1
end

3.相关结论的可视化

这里介绍可视化函数plotDecisionBoundary(theta, X, y);的编程思路，这个函数的定义如下：

function plotDecisionBoundary(theta, X, y)
plotData(X(:,2:3), y);%这个函数就是按类绘制散点图，其定义及原理见上文
hold on %锁定画布
if size(X, 2) <= 3 %这个判断使得这个函数变成通用函数，只有两个变量时画直线
    plot_x = [min(X(:,2))-2,  max(X(:,2))+2];
    %横坐标为x1，这里选取两个点确定一条直线
    plot_y = (-1./theta(3)).*(theta(2).*plot_x + theta(1));
    %这个计算就是西格蒙函数=0.5的临界线，具体推理见笔记
    plot(plot_x, plot_y)   
    % 绘制直线
    legend('Admitted', 'Not admitted', 'Decision Boundary')
    %标注图例
    axis([30, 100, 30, 100])
    %规定坐标取值范围
else
	%大于两个变量时，绘制等高线来作为边界范围
    % Here is the grid range
    u = linspace(-1, 1.5, 50);
    %横坐标取值范围-1到1.5分50等份
    v = linspace(-1, 1.5, 50);
    %纵坐标如是
    z = zeros(length(u), length(v));
    %所谓z就是西格蒙函数的参数z，即theta*X，这里先初始化为0，再逐一计算z
    for i = 1:length(u)
        for j = 1:length(v)
            z(i,j) = mapFeature(u(i), v(j))*theta;
            %计算每一个z，注意，这里的theta是多项式展开的，所以这里也要展开
            %测试的时候也要展开
        end
    end
    z = z'; % important to transpose z before calling contour
    %绘制等高线时要转置
    contour(u, v, z, [0, 0], 'LineWidth', 2)
    %第四个参数为level，意为在何等级处划线
    %如果要在某一特定层级k上划线，要定义为二元向量[k,k]
    %这里要在西格蒙函数值为0.5时，也就是z=0时划线，所以是[0,0]
end
hold off %解锁
end

这样，就完成了可视化
除了将一些变量进行可视化之外，这里再介绍一种评价算法好坏的评价函数
说白了就是将每个模型估计出来的输出与label相减，然后除以样本容量再*100
因为Logistic回归输出只有0和1所以这种方法尤其适用
评价方法代码如下：

mean(double(p == y)) * 100)
%逻辑判断转double数字，如果p==y则相当于1
%否则是0

4.多项式Logistic回归

如果图像的数据较为复杂不能简单地认为是画一条直线就能分类，就要考虑他们的多项式相乘的因子带来的影响，首先要将两个变量相乘得到多项式，再将多项式的每一项当做一个一次变量，这个就形成了多元一次Logistic回归

为了实现上述方法，定义函数如下：

function out = mapFeature(X1, X2)
degree = 6;%自由度，也就是最高次数
out = ones(size(X1(:,1)));
%定义一个用于储存多元变量，并且顺便把x0填上便于后面的计算
for i = 1:degree
    for j = 0:i
        out(:, end+1) = (X1.^(i-j)).*(X2.^j);
        添加多元变量
    end
end
end

这样就完成了对变量的扩展，有了这样的函数就可以对数据进行预处理了
但是要注意的是，多项式回归与刚才的回归不同的地方就是要进行正则化，所以代价函数和梯度会有所不同，后面我们详述
初始化过程这里不再赘述

data = load('ex2data2.txt');
X = data(:, [1, 2]); y = data(:, 3);
plotData(X, y);
hold on;
xlabel('Microchip Test 1')
ylabel('Microchip Test 2')
legend('y = 1', 'y = 0')
hold off;
X = mapFeature(X(:,1), X(:,2));
initial_theta = zeros(size(X, 2), 1);
lambda = 1;

计算代价函数和梯度

这里要注意，正则化中x0是不参与的，也就要在弄一个不含x0的X与含x0的X混合使用进行计算

function [J, grad] = costFunctionReg(theta, X, y, lambda)
m = length(y); % 样本容量
theta_1=[0;theta(2:end)];    % 弄一个不含x0的，待会x0不参与的部分用这个
J= -1 * sum( y .* log( sigmoid(X*theta) ) + (1 - y ) ...
	.* log( (1 - sigmoid(X*theta)) ) ) / m  + lambda/(2*m) * theta_1' * theta_1; 
grad = ( X' * (sigmoid(X*theta) - y ) )/ m + lambda/m * theta_1; 
end

剩下的就与一开始的如出一辙，现演示如下：

X = mapFeature(X(:,1), X(:,2));
[cost, grad] = costFunctionReg(initial_theta, X, y, lambda);
initial_theta = zeros(size(X, 2), 1);
options = optimset('GradObj', 'on', 'MaxIter', 400);
[theta, J, exit_flag] = ...
	fminunc(@(t)(costFunctionReg(t, X, y, lambda)), initial_theta, options);
plotDecisionBoundary(theta, X, y);
hold on;
title(sprintf('lambda = %g', lambda))
xlabel('Microchip Test 1')
ylabel('Microchip Test 2')
legend('y = 1', 'y = 0', 'Decision boundary')
hold off;
p = predict(theta, X);

调整λ还能演示欠拟合和过拟合等情况
这是λ=0的情况

这是λ=100的情况