[吴恩达机器学习exercise1：线性回归]

本文链接：https://blog.youkuaiyun.com/haronchou/article/details/104991251

1、单变量线性回归

根据房屋的size预测房屋的价格。

data = load('ex1data1.txt');%txt格式的数据，可以读入到Mat矩阵中！！
X = data(:, 1); y = data(:, 2);
m = length(y); % number of training examples

plot(X, y, 'rx', 'MarkerSize', 10); % Plot the data
ylabel('Profit in $10,000s'); % Set the y   axis label
xlabel('Population of City in 10,000s'); % Set the x   axis label

$\small h_{\theta}(x)=\theta_0+\theta_1x_1$ ，只有一个变量，用一条直线去拟合数据，得到参数θ0，θ1
如何评价这条直线拟合的好坏呢？

代价函数/拟合误差/建模误差

最小均方误差/最小二乘误差是回归问题常用的手段，使得所有样本的建模误差的平方和最小的模型参数

$\small J(\theta_0,\theta_1)$ 是关于θ0和θ1的凸函数，目标是找到合适的θ0和θ1使得建模误差/代价函数J最小，那么我们的模型 $\small h_{\theta}(x)=\theta_0+\theta_1x_1$ 就是最好的。就找到了最优的直线进行拟合。

梯度下降求解上述的代价函数，得到最优解θ0和θ1
- 开始随机选择一个参数组合（θ0, θ1, ... , θn）计算代价函数
- 选择一个让代价函数J下降最多的参数组合，持续到一个局部最小值
问题：(a)没有尝试完所有的参数组合，可能陷入局部最小值而非全局最小值 (b) 选择不同的初始参数组合，可能会找到不同的局部最小值
批量梯度下降 batch gradient descent：
- 在梯度下降的每一步中，都用到了所有的样本
- 求导时，要对所有m个训练样本求和
- α学习率——决定lost function沿着下降程度最大的方向迈进的步长有多大

也有不考虑整个训练集，只关注一些样本的小集合。此时就不是批量梯度下降，如随机梯度下降。

同时更新：vs 交替更新
同时更新需要一个临时变量，这样才能保证对θ0和θ1的同步更新，否则就是交替更新。

学习率α的选择
太小：下降慢，需要多步迭代才能达到全局最低点
太大：可能跳过最低点，无法收敛，甚至发散
但即使α固定，梯度下降也能收敛到局部最小，为什么？

$\small \theta_1=\theta_1-\alpha \frac{\mathrm{d} }{\mathrm{d} \theta_1}J(\theta_1)$ 越接近local minimum时，梯度 $\small \frac{\mathrm{d} }{\mathrm{d} \theta_1}J(\theta_1)$ 的值在慢慢变小，自然每次的步长都在减小。

在程序中， $\small X=\begin{bmatrix} X_1\\ X_2 \end{bmatrix}$ ，每一行为一个样本，每一列为特征。m为样本个数。

对于一个样本而言，其假设函数为 ${\color{Red} }h_\theta(x)=\theta_0+\theta_1x_1+...+\theta_nx_n=\theta^Tx$

矢量化表示代价函数为： $\small J=\frac{1}{2m}\sum(X*\theta-y)^2$ ，每行为一个样本的代价函数，对列求和，得到所有样本的代价函数之和。同时梯度下降的矢量化表示为： $\small \theta_j=\theta_j-\alpha\sum_{i=1}^{m}(X*\theta-y).*X(:,j)$ 。

X = [ones(m, 1), data(:,1)]; % Add a column of ones to x
theta = zeros(2, 1); % initialize fitting parameters

% Some gradient descent settings
iterations = 1500;
alpha = 0.01;

m = length(y); % number of training examples
J = sum((X * theta - y).^2) / (2*m);     % X(79,2)  theta(2,1)

%% 计算梯度下降
J_history = zeros(iterations , 1);
theta_s = theta;%初始theta值为0，从0开始梯度下降

for iter = 1:iterations % 固定迭代次数，假定在此之内一定收敛
    % 注意矢量化的方式
    theta(1) = theta(1) - alpha / m * sum(X * theta_s - y);       
    theta(2) = theta(2) - alpha / m * sum((X * theta_s - y) .* X(:,2));    
     % 必须同时更新theta(1)和theta(2)，所以不能用X * theta,而要用theta_s存储上次结果。
    theta_s = theta; 
    J_history(iter) = sum((X * theta - y).^2) / (2*m);
end

% Plot the linear fit
hold on; % keep previous plot visible
plot(X(:,2), X*theta, '-')
legend('Training data', 'Linear regression')
hold off 

predict1 = [1, 3.5] *theta;
predict2 = [1, 7] * theta;

Theta found by gradient descent: -3.630291 1.166362

For population = 35,000, we predict a profit of 4519.767868
For population = 70,000, we predict a profit of 45342.450129

可视化代价函数J：选择不同的θ0和θ1，得到不同的J。surf用于查看三维图，contour用于画等高线图，并标记出最优解的位置。

% Grid over which we will calculate J
theta0_vals = linspace(-10, 10, 100);
theta1_vals = linspace(-1, 4, 100);

% initialize J_vals to a matrix of 0's
J_vals = zeros(length(theta0_vals), length(theta1_vals));

% Fill out J_vals
for i = 1:length(theta0_vals)
    for j = 1:length(theta1_vals)
	  t = [theta0_vals(i); theta1_vals(j)];    
	  J_vals(i,j) = sum((X * t- y).^2) / (2*m); ;
    end
end

% 由于meshgrids工作在surf命令下，在surf之前需要转置J
J_vals = J_vals';
% Surface plot
figure;surf(theta0_vals, theta1_vals, J_vals)
xlabel('\theta_0'); ylabel('\theta_1');

% Contour plot
figure;
% Plot J_vals as 15 contours spaced logarithmically between 0.01 and 100
contour(theta0_vals, theta1_vals, J_vals, logspace(-2, 3, 20))
xlabel('\theta_0'); ylabel('\theta_1');
hold on;
plot(theta(1), theta(2), 'rx', 'MarkerSize', 10, 'LineWidth', 2);

2、多变量线性回归

两个变量的回归，首先要进行feature scaling，即数据标准化处理。

多变量时的梯度下降： $\theta=\theta-\alpha/m{\color{Red} X'(X*\theta-y)}$ 。后面的部分十分重要，矢量化的表示方式。

theta的值为：

340412.659574
110631.050279
-6649.474271

%% Load Data
data = load('ex1data2.txt');
X = data(:, 1:2);
y = data(:, 3);
m = length(y);

% 特征标准化处理
mu = mean(X);       %  mean value 对列求均值
sigma = std(X);     %  standard deviation 对列求标准差
X= (X - repmat(mu,size(X,1),1)) ./  repmat(sigma,size(X,1),1);

% Add intercept term to X
X = [ones(m, 1) X];

% 梯度下降求解
% Choose some alpha value
alpha = 0.01;
num_iters = 8500;

% Init Theta and Run Gradient Descent 
theta = zeros(3, 1);
m = length(y); % number of training examples
J_history = zeros(num_iters, 1);

for iter = 1:num_iters
     theta = theta - alpha / m * X' * (X * theta - y); 
    J_history(iter) = J = sum((X * theta - y).^2) / (2*m);    
end

对新数据进行预测：

price = [1 (([1650 3]-mu) ./ sigma)] * theta ;

2.1 正规方程求解

同样的，把每一行作为一个样本，

$X=\begin{bmatrix} 1 & X_1^T\\ 1&X_2^T \\ ... & ...\\ 1&X_m^T \end{bmatrix}$

那么代价函数可以写为： $J=\frac{1}{2m}\sum_{i=1}^{m}(X*\theta -y)^2=\frac{1}{2m}(X*\theta-y)^T(X*\theta-y)$ ，代价函数对θ求导的结果为：

$\frac{\partial J}{\partial \theta }=X^T(X\Theta-y)$ ，当这个偏导数为0时，取得最优解为： $\hat{\theta}=(X^TX)^{-1}X^Ty$

前提是可逆，即X'X要可逆，才能用正规方程的方式求解。X'X就是相关函数，相关函数是实对称阵，当X的特征是不相关联的时候，X‘X是实对称阵，一定可以相似对角化，一定是非奇异帧，一定存在逆。

theta = pinv( X' * X ) * X' * y;
price = [1 1650 3] * theta ;