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

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:\small \theta_j=\theta_j-\alpha\sum_{i=1}^{m}(X*\theta-y).*X(:,j)
    • 在梯度下降的每一步中,都用到了所有的样本
    • 求导时,要对所有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)}。后面的部分十分重要,矢量化的表示方式。

代价函数J的收敛情况

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 ;

 

  • 0
    点赞
  • 4
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值