2、逻辑回归

【逻辑回归函数模型】
1、训练数据绘图
cd D:\study\AI\data\ex2
data = load('ex2data1.txt');
X = data(:, [1, 2]); 
y = data(:, 3);

pos = find(y==1);  % y取1的所有行
plot(X(pos, 1), X(pos, 2), 'k+','LineWidth', 2, 'MarkerSize', 7);
hold on;
neg = find(y == 0);  % y取0的所有行
plot(X(neg, 1), X(neg, 2), 'ko', 'MarkerFaceColor', 'y', 'MarkerSize', 7);
xlabel('Exam 1 score'); 
ylabel('Exam 2 score'); 
% Specified in plot order
legend('Admitted', 'Not admitted'); 

2、计算损失函数和梯度
[m, n] = size(X);
X = [ones(m, 1) X];
initial_theta = zeros(n + 1, 1);
test_theta = [-24; 0.2; 0.2];
[cost, grad] = costFunction(test_theta, X, y);
fprintf('\nCost at test theta: %f\n', cost);
fprintf('Gradient at test theta: \n');
fprintf(' %f \n', grad);

3、高级算法
%  Set options for fminunc
options = optimset('GradObj', 'on', 'MaxIter', 400);
%  This function will return theta and the cost 
[theta, cost] = fminunc(@(t)(costFunction(t, X, y)), initial_theta, options);
% Print theta to screen
fprintf('Cost at theta found by fminunc: %f\n', cost);
fprintf('Expected cost (approx): 0.203\n');
fprintf('theta: \n');
fprintf(' %f \n', theta);
fprintf('Expected theta (approx):\n');
fprintf(' -25.161\n 0.206\n 0.201\n');

4、预测
prob = sigmoid([1 45 85] * theta);
fprintf(['For a student with scores 45 and 85, 'probability of %f\n'], prob);
fprintf('Expected value: 0.775 +/- 0.002\n\n');
p = predict(theta, X);
fprintf('Train Accuracy: %f\n', mean(double(p == y)) * 100);
fprintf('Expected accuracy (approx): 89.0\n');
fprintf('\n');

【正则化】
1、训练数据绘图
cd D:\study\AI\data\ex2
data = load('ex2data2.txt');
X = data(:, [1, 2]); 
y = data(:, 3);

pos = find(y==1);  % y取1的所有行
plot(X(pos, 1), X(pos, 2), 'k+','LineWidth', 2, 'MarkerSize', 7);
hold on;
neg = find(y == 0);  % y取0的所有行
plot(X(neg, 1), X(neg, 2), 'ko', 'MarkerFaceColor', 'y', 'MarkerSize', 7);
xlabel('Exam 1 score'); 
ylabel('Exam 2 score'); 
% Specified in plot order
legend('Admitted', 'Not admitted'); 

2、正则化计算
% 多项式特征映射
X = mapFeature(X(:,1), X(:,2));
initial_theta = zeros(size(X, 2), 1);
[cost, grad] = costFunctionReg(initial_theta, X, y, 1);
fprintf('Cost at initial theta (zeros): %f\n', cost);
fprintf(' %f \n', grad);

test_theta = ones(size(X,2),1);
[cost, grad] = costFunctionReg(test_theta, X, y, 10);
fprintf('\nCost at test theta (with lambda = 10): %f\n', cost);
fprintf(' %f \n', grad);

3、高级算法实现正则化
lambdas = [0, 1, 10, 100];
for lambda = lambdas,
    test(X, y, lambda);
end;
function test = test(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);
fprintf('Train Accuracy: %f\n', mean(double(p == y)) * 100);
fprintf('Expected accuracy (with lambda = 1): 83.1 (approx)\n');
end

【工具函数】
% logistic函数
function g = sigmoid(z)
g = zeros(size(z));
g = 1./(1 + exp(-z));
end;
% 损失函数和梯度
function [J, grad] = costFunction(theta, X, y)
m = length(y);
J = 0;
grad = zeros(size(theta));
% J(θ) = −y * log(g(x)) − (1 − y) * log(1 − g(x))
J = (-y'*log(sigmoid(X*theta)) - (1-y)'*(log(1-sigmoid(X*theta))))/m;
% ▽J(θ)= x′*(g(θ)-y)/m
grad = X'*(sigmoid(X*theta) - y)/m;
end;
% 正则约束损失函数和梯度
function [J, grad] = costFunctionReg(theta, X, y, lambda)
m = length(y);
grad = zeros(size(theta));
J = 1/m * (-y' * log(sigmoid(X*theta)) - (1 - y') * log(1 - sigmoid(X * theta))) + lambda/2/m*sum(theta(2:end).^2);
grad(1,:) = 1/m * (X(:, 1)' * (sigmoid(X*theta) - y));
grad(2:size(theta), :) = 1/m * (X(:, 2:size(theta))' * (sigmoid(X*theta) - y)) + lambda/m*theta(2:size(theta), :);
end;
% 预测函数
function p = predict(theta, X)
m = size(X, 1); % Number of training examples
p = zeros(m, 1);
p = sigmoid(X * theta)>=0.5;
end;
% 多项式特征映射
function out = mapFeature(X1, X2)
degree = 6;
out = ones(size(X1(:,1)));
for i = 1:degree
    for j = 0:i
        out(:, end+1) = (X1.^(i-j)).*(X2.^j);
    end
end
end;
% 决策边界绘图
function plotDecisionBoundary(theta, X, y)
pos = find(y==1);  % y取1的所有行
plot(X(pos, 2), X(pos, 3), 'k+','LineWidth', 2, 'MarkerSize', 7);
hold on;
neg = find(y == 0);  % y取0的所有行
plot(X(neg, 2), X(neg, 3), 'ko', 'MarkerFaceColor', 'y', 'MarkerSize', 7);
xlabel('Exam 1 score'); 
ylabel('Exam 2 score'); 
% Specified in plot order
legend('Admitted', 'Not admitted'); 
hold on
% 线性一阶决策边界直接绘制曲线图
if size(X, 2) <= 3
    plot_x = [min(X(:,2))-2,  max(X(:,2))+2];
    plot_y = (-1./theta(3)).*(theta(2).*plot_x + theta(1));
    plot(plot_x, plot_y)
    legend('Admitted', 'Not admitted', 'Decision Boundary')
    axis([30, 100, 30, 100])
else % 高阶决策边界通过等高线绘制Z=0切线图
    u = linspace(-1, 1.5, 50);
    v = linspace(-1, 1.5, 50);
    z = zeros(length(u), length(v));
    for i = 1:length(u)
        for j = 1:length(v)
            % 多项式特征映射的结构和X完全一样
            z(i,j) = mapFeature(u(i), v(j))*theta;
        end
    end
    z = z';
    contour(u, v, z, [0, 0], 'LineWidth', 2)
end
hold off
end;