正则化的线性回归以及偏差VS方差

正则化的线性回归

背景：数据集中包含水位变化的历史记录x和水坝的水量y。

任务一 可视化数据集

我们将数据集分为三部分：

• 训练集：X，y
• 交叉验证集：Xval，yval
• 测试集：Xtest，ytest

因此，本任务只需将训练集可视化即可。在ex5.m文件已将该任务代码准备好了，我们只需运行即可：

% Load from ex5data1: 
% You will have X, y, Xval, yval, Xtest, ytest in your environment
load ('ex5data1.mat');

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

% Plot training data
plot(X, y, 'rx', 'MarkerSize', 10, 'LineWidth', 1.5);
xlabel('Change in water level (x)');
ylabel('Water flowing out of the dam (y)');

其运行结果为：

任务二 线性回归的正则化代价函数

我们先将正则化代价函数公式列出：

其向量化后公式为：
　J(θ) = ((Xθ - y)^T(Xθ - y) + λθ_t^Tθ_t) / 2m
其中，θ_t表示将θ的第一列替换为0，即θ₀不参与正则化操作。

因此，我们可根据上述公式在linearRegCostFunction.m文件中键入如下代码：

theta_1 = [0; theta(2:end)];
J = ((X * theta - y)' * (X * theta - y)) / (2 * m) + lambda / (2 * m) * theta_1' * theta_1;

任务三 线性回归的正则化下降梯度

同样的，我们先将公式列出：

其向量化后的公式为：
　grad = (X^T(Xθ - y) + λθ_t) / m
其中θ_t同上。

因此，在linearRegCostFunction.m文件中继续键入如下代码：

grad = (X' * (X * theta - y) + lambda * theta_1) / m;

任务四 拟合线性回归

一旦你的代价函数和下降梯度运行正常，下一步就是在ex5.m文件调用并运行trainLinearReg.m文件中的代码，通过使用fmincg函数计算出使得代价函数最小化的θ。

在该任务中参数θ为2维向量，因此我们将正则化参数λ的值设为0。为何将正则化参数λ的值设为0？这是因为正则化对于低维度的θ没有太大的帮助。

trainLinearReg.m文件中的代码如下：

function [theta] = trainLinearReg(X, y, lambda)
%TRAINLINEARREG Trains linear regression given a dataset (X, y) and a
%regularization parameter lambda
%   [theta] = TRAINLINEARREG (X, y, lambda) trains linear regression using
%   the dataset (X, y) and regularization parameter lambda. Returns the
%   trained parameters theta.
%

% Initialize Theta
initial_theta = zeros(size(X, 2), 1); 

% Create "short hand" for the cost function to be minimized
costFunction = @(t) linearRegCostFunction(X, y, t, lambda);

% Now, costFunction is a function that takes in only one argument
options = optimset('MaxIter', 200, 'GradObj', 'on');

% Minimize using fmincg
theta = fmincg(costFunction, initial_theta, options);

end

ex5.m文件中该部分代码如下：

%  Train linear regression with lambda = 0
lambda = 0;
[theta] = trainLinearReg([ones(m, 1) X], y, lambda);

%  Plot fit over the data
plot(X, y, 'rx', 'MarkerSize', 10, 'LineWidth', 1.5);
xlabel('Change in water level (x)');
ylabel('Water flowing out of the dam (y)');
hold on;
plot(X, [ones(m, 1) X]*theta, '--', 'LineWidth', 2)
hold off;

fprintf('Program paused. Press enter to continue.\n');
pause;

运行结果为：

偏差与方差

高偏差的模型通常对训练集的拟合不太好，即欠拟合问题；高方差的模型通常对训练集的拟合非常完美，但对于交叉验证集或测试集的拟合不太好，即过拟合问题。

因此，该小节将练习绘制学习曲线来诊断偏差与方差的问题。

任务一 学习曲线

为了绘制学习曲线，我们需要计算出J_train(θ)和J_CV(θ)。

其中J_train(θ)的计算公式为：

J_CV(θ)的计算公式为：

因此，我们先需要利用trainLinearReg函数计算出使得代价函数最下化的θ的值；然后在使用linearRegCostFunction函数分别计算J_train(θ)和J_CV(θ)。

注：在使用linearRegCostFunction函数时，要注意将正则化参数λ = 0。

learningCurve.m文件中的具体代码如下：

for i = 1 : m
    theta = trainLinearReg(X(1:i, :), y(1:i), lambda);
    error_train(i) = linearRegCostFunction(X(1:i, :), y(1:i), theta, 0);
    error_val(i) = linearRegCostFunction(Xval, yval, theta, 0);
end

该部分的运行结果为：

多项式回归

在之前的部分，我们的线性模型对数据的拟合不太好，即出现欠拟合问题。在本小节，我们通过增加特征变量来解决欠拟合问题。

对于多项式回归，我们的假设函数h_θ(x)为：
h_θ(x) = θ₀ + θ₁ * (waterLevel) + θ₂ * (waterLevel)² + ... + θ_p * (waterLevel)^p

现在，我们需要在数据集中增添高阶幂的特征变量。因此，我们需要在polyFeatures.m文件中键入相关代码，使得数据集X变为一个m*p的矩阵。

polyFeatures.m文件的相关代码如下：

for i = 1 : p
    X_poly(:, i) = X .^ i;
end

任务一 学习多项式回归

对于该部分的练习，我们使用8次幂的多项式回归模型。由于特征变量在多项式回归模型中，其取值范围各不相同。因此，我们需要对特征变量归一化。

featureNormalize.m文件中特征变量归一化代码如下：

function [X_norm, mu, sigma] = featureNormalize(X)
%FEATURENORMALIZE Normalizes the features in X 
%   FEATURENORMALIZE(X) returns a normalized version of X where
%   the mean value of each feature is 0 and the standard deviation
%   is 1. This is often a good preprocessing step to do when
%   working with learning algorithms.

mu = mean(X);
X_norm = bsxfun(@minus, X, mu);

sigma = std(X_norm);
X_norm = bsxfun(@rdivide, X_norm, sigma);


% ============================================================

end

然后，我们将正则化参数λ = 0，利用trainLinearReg函数计算出使得代价函数最下化的θ的值。

最后，我们利用linearRegCostFunction函数分别计算Jtrain(θ)和JCV(θ)，绘制出学习曲线。

该部分代码如下：

%% =========== Part 6: Feature Mapping for Polynomial Regression =============
%  One solution to this is to use polynomial regression. You should now
%  complete polyFeatures to map each example into its powers
%

p = 8;

% Map X onto Polynomial Features and Normalize
X_poly = polyFeatures(X, p);
[X_poly, mu, sigma] = featureNormalize(X_poly);  % Normalize
X_poly = [ones(m, 1), X_poly];                   % Add Ones

% Map X_poly_test and normalize (using mu and sigma)
X_poly_test = polyFeatures(Xtest, p);
X_poly_test = bsxfun(@minus, X_poly_test, mu);
X_poly_test = bsxfun(@rdivide, X_poly_test, sigma);
X_poly_test = [ones(size(X_poly_test, 1), 1), X_poly_test];         % Add Ones

% Map X_poly_val and normalize (using mu and sigma)
X_poly_val = polyFeatures(Xval, p);
X_poly_val = bsxfun(@minus, X_poly_val, mu);
X_poly_val = bsxfun(@rdivide, X_poly_val, sigma);
X_poly_val = [ones(size(X_poly_val, 1), 1), X_poly_val];           % Add Ones

fprintf('Normalized Training Example 1:\n');
fprintf('  %f  \n', X_poly(1, :));

fprintf('\nProgram paused. Press enter to continue.\n');
pause;



%% =========== Part 7: Learning Curve for Polynomial Regression =============
%  Now, you will get to experiment with polynomial regression with multiple
%  values of lambda. The code below runs polynomial regression with 
%  lambda = 0. You should try running the code with different values of
%  lambda to see how the fit and learning curve change.
%

lambda = 0;
[theta] = trainLinearReg(X_poly, y, lambda);

% Plot training data and fit
figure(1);
plot(X, y, 'rx', 'MarkerSize', 10, 'LineWidth', 1.5);
plotFit(min(X), max(X), mu, sigma, theta, p);
xlabel('Change in water level (x)');
ylabel('Water flowing out of the dam (y)');
title (sprintf('Polynomial Regression Fit (lambda = %f)', lambda));

figure(2);
[error_train, error_val] = ...
    learningCurve(X_poly, y, X_poly_val, yval, lambda);
plot(1:m, error_train, 1:m, error_val);

title(sprintf('Polynomial Regression Learning Curve (lambda = %f)', lambda));
xlabel('Number of training examples')
ylabel('Error')
axis([0 13 0 100])
legend('Train', 'Cross Validation')

fprintf('Polynomial Regression (lambda = %f)\n\n', lambda);
fprintf('# Training Examples\tTrain Error\tCross Validation Error\n');
for i = 1:m
    fprintf('  \t%d\t\t%f\t%f\n', i, error_train(i), error_val(i));
end

fprintf('Program paused. Press enter to continue.\n');
pause;

运行结果为：

学习曲线

任务二 调整正则化参数（选做）

当正则化参数λ = 1时，其运行结果为：

学习曲线

当正则化参数λ = 100时，其运行结果为：

学习曲线

任务三 通过交叉验证集选择正则化参数

当正则化参数λ∈{0, 0.001, 0.003, 0.03, 0.1, 0.3, 1, 3, 10}时，分别计算出J_train(θ)和J_CV(θ)。

在validationCurve.m文件中键入如下代码：

for i = 1 : length(lambda_vec)
    lambda = lambda_vec(i);
    theta = trainLinearReg(X, y, lambda);
    error_train(i) = linearRegCostFunction(X, y, theta, 0);
    error_val(i) = linearRegCostFunction(Xval, yval, theta, 0);
end

然后，ex5.m文件中的相关代码通过J_train(θ)和J_CV(θ)的值，绘制出相关函数图。

[lambda_vec, error_train, error_val] = ...
    validationCurve(X_poly, y, X_poly_val, yval);

close all;
plot(lambda_vec, error_train, lambda_vec, error_val);
legend('Train', 'Cross Validation');
xlabel('lambda');
ylabel('Error');

fprintf('lambda\t\tTrain Error\tValidation Error\n');
for i = 1:length(lambda_vec)
    fprintf(' %f\t%f\t%f\n', ...
            lambda_vec(i), error_train(i), error_val(i));
end

fprintf('Program paused. Press enter to continue.\n');
pause;

其运行结果为：

任务四 计算测试集误差（选做）

在实际开发中，除了计算J_train(θ)和J_CV(θ)，我们还需计算J_test(θ)。

参考代码：

for i = 1 : m
    theta = trainLinearReg(X(1:i, :), y(1:i), lambda);
    error_test(i) = linearRegCostFunction(Xtest(1:i, :), ytest(1:i), theta, 0);
end

任务五 绘制随机数据集的学习曲线（选做）

从数据集中随机抽取60%的数据作为训练集，20%的数据作为交叉验证集和20%的数据作为测试集。由于此部分答案不唯一就不过多叙述。