线性回归、梯度下降算法 线性回归、梯度下降算法

中介 ——> 模型
从业 ——> 学习
影响 ——> 权重
房子 ——> 样本
价格 ——> 输出

fprintf('Plotting Data ...
')
data = load('ex1data1.txt');
X = data(:, 1); 
y = data(:, 2);
m = length(y); % number of training examples

% Plot Data
% Note: You have to complete the code in plotData.m
plotData(X, y);

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

function J = computeCost(X, y, theta)
m = length(y); % number of training examples

% J=1/(2*m)*sum((X*theta-y).^2); 

J=1/(2*m)*(X*theta-y)'*(X*theta-y);

end

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

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

fprintf('
Testing the cost function ...
')
% compute and display initial cost
J = computeCost(X, y, theta);
fprintf('With theta = [0 ; 0]
Cost computed = %f
', J);
fprintf('Expected cost value (approx) 32.07
');

% further testing of the cost function
J = computeCost(X, y, [-1 ; 2]);
fprintf('
With theta = [-1 ; 2]
Cost computed = %f
', J);
fprintf('Expected cost value (approx) 54.24
');

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

function [theta, J_history] = gradientDescent(X, y, theta, alpha, num_iters)
m = length(y); % number of training examples
J_history = zeros(num_iters, 1);  % one column

for iter = 1:num_iters
    grad = 1/m * X'*(X*theta-y);
    theta = theta - alpha * grad; 
    J_history(iter) = computeCost(X, y, theta);
end
endfprintf('
Running Gradient Descent ...
')
% run gradient descent
theta = gradientDescent(X, y, theta, alpha, iterations);

% print theta to screen
fprintf('Theta found by gradient descent:
');
fprintf('%f
', theta);
fprintf('Expected theta values (approx)
');
fprintf(' -3.6303
  1.1664

');

fprintf('
Running Gradient Descent ...
')
% run gradient descent
theta = gradientDescent(X, y, theta, alpha, iterations);

% print theta to screen
fprintf('Theta found by gradient descent:
');
fprintf('%f
', theta);
fprintf('Expected theta values (approx)
');
fprintf(' -3.6303
  1.1664

');

% Predict values for population sizes of 35,000 and 70,000

predict1 = [1, 3.5] *theta;
fprintf('For population = 35,000, we predict a profit of %f
',...
    predict1*10000);
    
predict2 = [1, 7] * theta;
fprintf('For population = 70,000, we predict a profit of %f
',...
    predict2*10000);

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

% Print out some data points
fprintf('First 10 examples from the dataset: 
');
fprintf(' x = [%.0f %.0f], y = %.0f 
', [X(1:10,:) y(1:10,:)]');

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

% Scale features and set them to zero mean
fprintf('Normalizing Features ...
');

% 将特征的值，都约束在[-1,1]区间内

[X mu sigma] = featureNormalize(X);

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

function [X_norm, mu, sigma] = featureNormalize(X)
X_norm = X;
mu = zeros(1, size(X, 2));   % get the count of column
sigma = zeros(1, size(X, 2));
[X_norm,mu,sigma] = zscore(X);
end

function [theta] = normalEqn(X, y)

theta = zeros(size(X, 2), 1);

theta = pinv(X'*X)*X'*y;

end

data = csvread('ex1data2.txt');
X = data(:, 1:2);
y = data(:, 3);
m = length(y);

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

% Calculate the parameters from the normal equation
theta = normalEqn(X, y);

% Display normal equation's result
fprintf('Theta computed from the normal equations: 
');
fprintf(' %f 
', theta);
fprintf('
');


% Estimate the price of a 1650 sq-ft, 3 br house
% ====================== YOUR CODE HERE ======================
% price = 0; % You should change this

price = 1650; % You should change this
price = [1,1650,3]*theta;

一、预测

先来看看这样一个场景：

假如你手头有一套房子要出售，你咨询了房产中介。中介跟你要了一系列的数据，例如房子面积、位置、楼层、年限等，然后进行一系列计算后，给出了建议的定价。

房产中介是如何帮你定价的？

“中介”通过他多年的“从业”经验，知道哪些因素会影响房子的价格，且知道各自的“影响”有多大，于是在接过“你的房子”时，他就能通过自已的经验计算出“价格”了。

当然，这个价格，不同的中介，得到的也不同。有经验的中介，可以很准确的预测到价格。而经验不足者，可以会有很大的偏差。

注意上面引号内的东西。将引号内的东西做如下抽象：

中介 ——> 模型
从业 ——> 学习
影响 ——> 权重
房子 ——> 样本
价格 ——> 输出

这便可以当作是一个线性回归模型。

二、线性回归

线性回归，用于解决数值预测问题。

1、模型假设

假设有一类数据样本y 表示该样本的值。

假设 x 之间，存在某种线性关系，那么数学上，可以用如下方程表示：

h(x)=θ0x0+θ1x1+...+θnxn

用矩阵表示为：

hθ(x)=θTX

显然，要确定 θ 的值便可以。

那么，θ 应该怎么求？

2、梯度下降算法

2.1 直观理解梯度下降

上述的问题，归结到求 θ。

求取 θ。

显然是OK的，但关键是 θ 每次要变化多大的量合适？当然了，每次就修改一点点，例如0.00000000001，然后修改了100年，终于得到一个合适的参数，这样理论上也是可以的。

那有没有什么办法，可以让 θ “快速的变化”到理想状态呢？

注意引号中的内容——快速的变化，想到什么没？是的，变化率！

数学上来讲，是求导！

用谁对谁求？当然是计算值与y之前误差了。

有了变化率，θ 要变大，还是变小？

梯度下降算法认为，往梯度反方向变化，即是做减法。

2.2 梯度下降法计算过程

误差函数

令误差函数如下式：

J=12m∑i=1m(hθ(x(i))−y(i))2

这个式子，解释一下。

误差用输出值与 y 值作差，这个没得说。平方是让差值恒为正。而前方的 1/2，是为了求导时将平方约去，简化运算。

求导

θ偏导数为

grad=1m∑i=1m(hθ(x(i))−y(i))xj(i)

θ

由梯度下降法，迭代计算α为学习率

θj:=θj−α∗grad

三、单一特征的线性回归

1、数据可视化

fprintf('Plotting Data ...
')
data = load('ex1data1.txt');
X = data(:, 1); 
y = data(:, 2);
m = length(y); % number of training examples

% Plot Data
% Note: You have to complete the code in plotData.m
plotData(X, y);

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

2、误差函数及测试

2.1、误差函数

根据公式，编写J函数

function J = computeCost(X, y, theta)
m = length(y); % number of training examples

% J=1/(2*m)*sum((X*theta-y).^2); 

J=1/(2*m)*(X*theta-y)'*(X*theta-y);

end

2.2、测试误差函数

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

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

fprintf('
Testing the cost function ...
')
% compute and display initial cost
J = computeCost(X, y, theta);
fprintf('With theta = [0 ; 0]
Cost computed = %f
', J);
fprintf('Expected cost value (approx) 32.07
');

% further testing of the cost function
J = computeCost(X, y, [-1 ; 2]);
fprintf('
With theta = [-1 ; 2]
Cost computed = %f
', J);
fprintf('Expected cost value (approx) 54.24
');

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

3、梯度下降

3.1、梯度下降算法

根据梯度下降公式，实现算法

function [theta, J_history] = gradientDescent(X, y, theta, alpha, num_iters)
m = length(y); % number of training examples
J_history = zeros(num_iters, 1);  % one column

for iter = 1:num_iters
    grad = 1/m * X'*(X*theta-y);
    theta = theta - alpha * grad; 
    J_history(iter) = computeCost(X, y, theta);
end
endfprintf('
Running Gradient Descent ...
')
% run gradient descent
theta = gradientDescent(X, y, theta, alpha, iterations);

% print theta to screen
fprintf('Theta found by gradient descent:
');
fprintf('%f
', theta);
fprintf('Expected theta values (approx)
');
fprintf(' -3.6303
  1.1664

');

3.2、运行梯度下降

传入训练集，计算出theta矩阵

fprintf('
Running Gradient Descent ...
')
% run gradient descent
theta = gradientDescent(X, y, theta, alpha, iterations);

% print theta to screen
fprintf('Theta found by gradient descent:
');
fprintf('%f
', theta);
fprintf('Expected theta values (approx)
');
fprintf(' -3.6303
  1.1664

');

4、预测

根据计算出来的theta矩阵，预测新的样本对应的值。

% Predict values for population sizes of 35,000 and 70,000

predict1 = [1, 3.5] *theta;
fprintf('For population = 35,000, we predict a profit of %f
',...
    predict1*10000);
    
predict2 = [1, 7] * theta;
fprintf('For population = 70,000, we predict a profit of %f
',...
    predict2*10000);

四、多特征的线性回归

多特征的线性回归，算法与单一特征基本一样，唯一不同的仅仅是X的特征数量、theta的参数数量。

另外，多特征的样本，各个特征的数据范围可能存在数量级的差别，最后得出的误差函数的等高线可能是一个不规则的图形。而梯度下降，是沿着误差函数的等高线的法向方向进行，对于一个不规则的图形，其法向方向是多变的，这样一来，梯度下降可能每次迭代的方向都会有差别，进而形成震荡，迭代速度也会变慢。

为使算法更快的找到最优解，使用“均值归一化”的进行标准化处理，将样本各个特征的数据，约束在相同的范围内。这样处理后，如果是二维特征，其误差函数的等高线就是一个近圆，梯度下降就直指圆心进行，能更快找到最优解。

1、标准化特征

1.1、 加载数据，并标准化处理

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

% Print out some data points
fprintf('First 10 examples from the dataset: 
');
fprintf(' x = [%.0f %.0f], y = %.0f 
', [X(1:10,:) y(1:10,:)]');

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

% Scale features and set them to zero mean
fprintf('Normalizing Features ...
');

% 将特征的值，都约束在[-1,1]区间内

[X mu sigma] = featureNormalize(X);

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

1.2、标准化函数（使用均值归一化算法）

计算方法：X_out = (X_in - mu)/sigma，其中，mu为样本均值，sigma为样本标准差

function [X_norm, mu, sigma] = featureNormalize(X)
X_norm = X;
mu = zeros(1, size(X, 2));   % get the count of column
sigma = zeros(1, size(X, 2));
[X_norm,mu,sigma] = zscore(X);
end

2、多维特征样本的误差函数、梯度下降算法

与一维样本的公式完全一样。不再详述。

五、正规方程解法

求解theta，除了使用梯度下降法外，还可以直接使用公式，进行计算。

上面讲到，如果是多维数据，在进行梯度下降时，为使计算速度加快，需要将多维特征数据进行标准化处理。而若是使用正规方程解法，则不需要，在原有的样本进行计算即可。

正规方程解法，因为需要求解逆矩阵，在样本数量很大的情况下，其运算速度会非常慢，所以，梯度下降法与正规方程解法，根据数据量的大小，进行选择。

经验数值是，当数据量大于1,000,000时，一般采用梯度下降法。

1、方程解法

function [theta] = normalEqn(X, y)

theta = zeros(size(X, 2), 1);

theta = pinv(X'*X)*X'*y;

end

2、预测

可以看到，正规方程解法，与梯度下降法，得出来的结果很接近。

data = csvread('ex1data2.txt');
X = data(:, 1:2);
y = data(:, 3);
m = length(y);

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

% Calculate the parameters from the normal equation
theta = normalEqn(X, y);

% Display normal equation's result
fprintf('Theta computed from the normal equations: 
');
fprintf(' %f 
', theta);
fprintf('
');


% Estimate the price of a 1650 sq-ft, 3 br house
% ====================== YOUR CODE HERE ======================
% price = 0; % You should change this

price = 1650; % You should change this
price = [1,1650,3]*theta;