关于 Octave 实现 logistic 回归的一些基本操作

前言

最近在 bilibili 看了点吴恩达的机器学习，感觉甚是经典，于是试着自己实现了一个简单的计算 logistic 回归的程序。

数据构造

生成 $[0,1]$ 范围内的随机数矩阵

rand(n); % 生成 n * n 的随机数方阵
rand(n,m); % 生成 n 行 m 列的随机数矩阵

生成 $10^4$ 个 $[0,1]\times [0,1]$ 中的散点

细节：此处不写 global 会给后期带来巨大的麻烦。

global N = 10000; % 数据集元素个数
global X = rand(2,N); % 生成一万个散点，每列是一个点，用来做数据集
global Y = zeros(N,1); % Y 用来储存标号
for i = 1:N,
	if X(1,i)^2 + X(2,i)^2 <= 1, % 圆内的点标号为 1, 圆外的点标号为 0
		Y(i) = 1;
	else；
		Y(i) = 0; % 可以不写，因为 Y 的初始值为 0
	end;
end;

我们称距离源点距离小于等于 $1$ 的点成为阳性点，大于 $1$ 的点称为阴性点，由此可以将点集分为两类。具体实践中拟合的数据一般都是在生产实践过程中产生的，此处使用随机数只是为了演示 logistic 回归的实现方式。

根据 $Y$ 向量筛选对应的点，并绘制散点图

setY = find(Y); % 筛选圆内数据
setN = find(1 .- Y); % 筛选圆外数据

hold on; % 在同一张图表上绘制多组数据
plot(X(1, setY), X(2, setY), 'b.'); % 蓝色显示圆内数据散点图
plot(X(1, setN), X(2, setN), 'r.'); % 红色显示圆外数据散点图
hold off;

到现在为止，数据集已经构造并检验完成了，接下来就是学习算法的主要内容。

logistic 回归

首先需要定义线性回归的代价函数，以及其梯度的计算方法。

logistic 假设函数 hypo

function jVal = hypo(theta, x) % 根据 logistic 假设计算预测值 y
	jVal = 1/(1+exp(-theta' * x));
end;

代价函数 costFunction

% 此时要保证 X 是 3 * N 的矩阵 (第一行全是 1)
function [jVal, gradient] = costFunction(theta) % theta 有三个维度 theta(1) 表示常数项
    global N;
    global X;
    global Y; % 必须要在独立的语句中声明对全局变量的使用 "global N, X, Y;" 不行
    
	jVal = 0;
	for i = 1: N, % 代价函数
		jVal -= (Y(i) * log(hypo(theta, X(:, i))) + (1 - Y(i)) * log(1 - hypo(theta, X(:, i))))/N; % logistic 代价函数
	end;
	
	gradient = zeros(3, 1);
	for i = 1: N,
		gradient += (hypo(theta, X(:, i)) - Y(i)) .* X(:, i) ./ N; % 计算梯度函数
	end;
end;

调用 octave 自带的优化函数 fminunc

X = [ones(1, N) ; X]; % 为 X 向量添加第一行数据 1
options = optimset("GradObj", "on", "MaxIter", 100); % 开启梯度函数，最多迭代 100 次
initialTheta = [0.5; 0.5; 0.5]; % theta 初始值
[optTheta, functionVal, exitFlag] = fminunc(@costFunction, initialTheta, options); % 运行优化算法

输出计算得到的分类直线，并检验分类效果

将数据点和直线绘制在同一个坐标系下

LX = 0: 0.001: 1; % 横坐标取值范围
LY = (-optTheta(1) .- optTheta(2) * LX) / optTheta(3); % 构造纵坐标 (下文见解释)
plot(X(2, setY), X(3, setY), "b.");
hold on; % 在同一张图上画图
plot(X(2, setN), X(3, setN), "r.");
plot(LX, LY, "k", "LineWidth", 2); % 黑色粗线条绘制分类线
hold off;

解释：由于我们采用 logistic 回归，最后训练得到的函数为：

$$h_{\theta }(x)=\frac{1}{1+e^{-{\theta }^{T}x}}$$

当 ${\theta }^{T}x=0$时，$h_{\theta }(x)=\frac{1}{2}$，因此，直线 ${\theta }^{T}x=0$ 就是阴性区域和阳性区域的分界线。写成直角坐标系下的方程为：${\theta }_1+{\theta }_{2}x+{\theta }_{3}y=0$，解得：

$$y=\frac{-{\theta }_1-{\theta }_{2}x}{{\theta }_3}$$

计算分类的正确率

predictY = zeros(N,1); % 计算预测结果
for i = 1: N,
	if hypo(optTheta, X(:, i)) > 0.5,
		predictY(i) = 1;
	else;
		predictY(i) = 0; % 由于初始值为 0，这句也可以不写
	end;
end;
acc = N - sum(xor(predictY, Y)); % 计算正确分类的元素个数
printf("acc = %.2f%%\n", acc / N * 100); % 输出正确率

最终效果

经过了大概 $30s$ 的训练，得到了命令行输出：

acc = 96.12%

绘图板输出：

无论是从图像还是从准确率上来看，这次 logistic 回归都还是比较成功的。