Exercise 2：Logistic Regression—實現一個邏輯迴歸
問題描述：用邏輯迴歸根據學生的考試成績來判斷該學生是否可以入學。

這裡的訓練資料(training instance)是學生的兩次考試成績，以及TA是否能夠入學的決定（y=0表示成績不合格，不予錄取；y=1表示錄取）

因此，需要根據trainging set 訓練出一個classification model。然後，拿著這個classification model 來評估新學生能否入學。

訓練資料的成績樣例如下：第一列表示第一次考試成績，第二列表示第二次考試成績，第三列表示入學結果（0–不能入學，1–可以入學）

34.62365962451697, 78.0246928153624,  0
30.28671076822607, 43.89499752400101, 0
35.84740876993872, 72.90219802708364, 0
60.18259938620976, 86.30855209546826, 1
...               ...

訓練資料圖形表示 如下：橫座標是第一次考試的成績，縱座標是第二次考試的成績，右上角的 + 表示允許入學，圓圈表示不允許入學。（分數決定命運，太悲慘了！）

該訓練資料的圖形 可以通過Matlab plotData函式畫出來,它呼叫Matlab中的plot函式和find函式，Matlab程式碼實現如下：
0. Visualizing the data

%Initialization
clear all; close all; clc

data = csvread('ex2data1.txt');
X = data(:, [1, 2]); y = data(:, 3);
%X讀取的是第一列和第二列的所有資料
%y讀取的是第三列的所有資料

錄取的點用+表示，未錄取的點用圓圈o表示。plotData函式的框架如下：

% ==================== Part 1: Plotting ====================

fprintf(['Plotting data with + indicating (y = 1) examples and o '
 
 ...
         'indicating (y = 0) examples.\n']);
function plotData(X, y)
    figure; hold on;
%====================== YOUR CODE HERE ======================
    xlabel('Exam 1 score')
    ylabel('Exam 2 score')
    legend('Admitted', 'Not admitted')
    hold off;
end

載入完資料之後，執行以下程式碼（呼叫自定義的plotData函式），將圖形畫出來：

function plotData(X, y)
figure; hold on;

% ====================== YOUR CODE HERE ======================


pos = find(y == 1);
neg = find(y == 0);

plot(X(pos,1),X(pos,2),'k+','Linewidth',2,'MarkerSize',7);
plot(X(neg,1),X(neg,2),'ko','MarkerFaceColor','y','MarkerSize',7);

註釋：
1.pos = find(y == 1）
找到y等於1的，資料的行數
neg = find(y == 0);
找到y等於0的，資料的行數
X(pos,1)表示X中，第pos行，第一列資料
X(pos,2)表示X中，第pos行，第二列資料
下面類似，不在贅述
2.plot(X(pos,1),X(pos,2),’k+’,’Linewidth’,2,’MarkerSize’,7);
表示。畫出，橫座標等於X(pos,1)縱座標等於X(pos,2)的點。用+號表示，線寬為2，大小為7
下面類似

圖形畫出來之後，對訓練資料就有了一個大體的視覺化的認識了。接下來就要實現 模型了，這裡需要訓練一個邏輯迴歸模型。

1. sigmoid function
要求實現如下函式，實現時z是向量，要對每個元素都做計算。

function g = sigmoid(z)
    g = zeros(size(z));設定初值
    g = 1 ./ (1 + exp(-z));
    %‘點除’ 表示 1 除以矩陣(向量)中的每一個元素
end

1. Cost function and gradient
   先來看一下計算前的準備工作

%% ============ Part 2: Compute Cost and Gradient ============

[m, n] = size(X);

X = [ones(m, 1) X];%向量化的常規套路，增加一行X0

initial_theta = zeros(n + 1, 1);%重置初值

[cost, grad] = costFunction(initial_theta, X, y);

X是個矩陣，它每列表示一個特徵，每行表示一個例項。theta是n+1個元素的列向量。y是m個元素的結果向量。

代價函式公式

可以先用一個變數hx求出h(x)，利用1中的函式

hx = sigmoid(X * theta);

求和號中每一項對應一行，可以使用向量化的方法計算，最後用sum求和。

J = (-1/m) * sum(y .* log(hx) + (1-y) .* log(1 - hx));

梯度的計算公式


代價函式costFunction（）的定義如下

function [J, grad] = costFunction(theta, X, y)
    m = length(y); % number of training examples
    grad = zeros(size(theta));
    hx = sigmoid(X * theta);
    J = (-1/m) * sum(y .* log(hx) + (1-y) .* log(1 - hx));

    n = length(theta);
    for j = 1 : n
        grad(j) = sum((hx - y) .* X(:, j)) / m;
         % X 為 training set 中的 feature variables, y 為training instance(訓練樣本的結果)結果，相當於對j
    end
end

通過呼叫costfunction.m檔案中定義的coustFunction函式，從而執行梯度下降演算法找到使代價函式J(theta)最小化的 邏輯迴歸模型引數theta。呼叫costFunction函式的程式碼如下：
3. fminunc

%  Set options for fminunc
options = optimset('GradObj', 'on', 'MaxIter', 400);
 %  Run fminunc to obtain the optimal theta
 %  This function will return theta and the cost 
[theta, cost] = ...
    fminunc(@(t)(costFunction(t, X, y)), initial_theta, options);
 % Plot Boundary
plotDecisionBoundary(theta, X, y);

註釋：
1.fminunc中的@(t)(costFunction(t, X, y))是為了把函式costFunction轉換為只接受一個引數t的函式。

2.從上面程式碼的最後一行可以看出，我們是通過 fminunc 呼叫 costFunction函式，來求得 theta的，而不是自己使用 Gradient descent 在for 迴圈求導來計算 theta。for迴圈中求導計算theta，
既然已經通過Gradient descent演算法求得了theta，將theta代入到假設函式中，就得到了 logistic regression model.

呼叫plotDecisionBoundary函式畫出點和決策邊界，其中相關的程式碼如下（兩個引數的情況）

if size(X, 2) <= 3
    % Only need 2 points to define a line, so choose two endpoints
    plot_x = [min(X(:,2))-2,  max(X(:,2))+2];

    % Calculate the decision boundary line
    plot_y = (-1./theta(3)).*(theta(2).*plot_x + theta(1));

    % Plot, and adjust axes for better viewing
    plot(plot_x, plot_y)

    % Legend, specific for the exercise
    legend('Admitted', 'Not admitted', 'Decision Boundary')
    axis([30, 100, 30, 100])

？？？？？？？？？？以上邊界函式程式碼看不懂，先放棄，以後再看，先往下進行。

ex2的第104行是不需要的，因為在上面函式裡已經添加了圖示。在執行legend(‘Admitted’, ‘Not admitted’)反而會覆蓋掉分界線的圖示。
將theta代入到假設函式中，就得到了 logistic regression model，用圖形表示如下：

4. predict 模型的評估（Evaluating logistic regression）

那如何估計，求得的邏輯迴歸模型是好還是壞呢？預測效果怎麼樣？因此，就需要拿一組資料測試一下，測試程式碼如下：

%% ============== Part 4: Predict and Accuracies ==============

prob = sigmoid([1 45 85] * theta); %這是一組測試資料，第一次考試成績為45，第二次成績為85
fprintf(['For a student with scores 45 and 85, we predict an admission ' ...
         'probability of %f\n\n'], prob);

% Compute accuracy on our training set
p = predict(theta, X);% 呼叫predict函式測試模型

fprintf('Train Accuracy: %f\n', mean(double(p == y)) * 100);

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

模型的測試結果如下：

For a student with scores 45 and 85, we predict an admission probability of 0.774323

Train Accuracy: 89.000000

那predict函式是如何實現的呢？predict.m 如下：

function p = predict(theta, X)

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

p = zeros(m, 1);

% ====================== YOUR CODE HERE ======================
% Instructions: Complete the following code to make predictions using
%               your learned logistic regression parameters. 
%               You should set p to a vector of 0's and 1's
%
p = X*theta >= 0;

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

end

備註：double(y)表示將引數y轉為雙精度浮搜尋點型別.mean(A)表示求A的平均值
非常簡單，只有一行程式碼：p = X * theta >= 0，原理如下：

當h(x)>=0.5時，預測y==1，而h(x)>=0.5 等價於 z>=0

Regularization
為什麼需要正則化？正則化就是為了解決過擬合問題(overfitting problem)。那什麼又是過擬合問題呢？

一般而言，當模型的特徵(feature variables)非常多，而訓練的樣本數目(training set)又比較少的時候，訓練得到的假設函式(hypothesis function)能夠非常好地匹配training set中的資料，此時的代價函式幾乎為0。下圖中最右邊的那個模型 就是一個過擬合的模型。

**所謂過擬合，從圖形上看就是：假設函式曲線完美地通過中樣本中的每一個點。也許有人會說：這不正是最完美的模型嗎？它完美地匹配了traing set中的每一個樣本呀！

過擬合模型不好的原因是：儘管它能完美匹配traing set中的每一個樣本，但它不能很好地對未知的 (新樣本例項)input instance 進行預測呀！通俗地講，就是過擬合模型的預測能力差。

因此，正則化(regularization)就出馬了。

前面提到，正是因為 feature variable非常多，導致 hypothesis function 的冪次很高，hypothesis function變得很複雜(彎彎曲曲的)，從而通過穿過每一個樣本點(完美匹配每個樣本)。如果新增一個”正則化項”，減少 高冪次的特徵變數的影響，那 hypothesis function不就變得平滑了嗎？

正如前面提到，梯度下降演算法的目標是最小化cost function，而現在把 theta(3) 和 theta(4)的係數設定為1000，設得很大，求偏導數時，相應地得到的theta(3) 和 theta(4) 就都約等於0了。

更一般地，我們對每一個theta(j)，j>=1，進行正則化，就得到了一個如下的代價函式：其中的 lambda(λ)就稱為正則化引數(regularization parameter)

從上面的J(theta)可以看出：如果lambda(λ)=0，則表示沒有使用正則化；如果lambda(λ)過大，使得模型的各個引數都變得很小，導致h(x)=theta(0)，從而造成欠擬合；如果lambda(λ)很小，則未充分起到正則化的效果。因此，lambda(λ)的值要合適。

最後，我們來看一個實際的過擬合的示例，原始的訓練資料如下圖：

lambda(λ)==1時，訓練出來的模型（hypothesis function）如下：Train Accuracy: 83.050847

lambda(λ)==0時，不使用正則化，訓練出來的模型（hypothesis function）如下：Train Accuracy: 87.288136

lambda(λ)==100時，訓練出來的模型（hypothesis function）如下：Train Accuracy: 61.016949

1. Visualizing the data
   這部分使用的程式碼和上面的相同
   
   但這次的資料不能用直線來擬合出邊界。
2. Feature mapping

使用已有的特徵的冪來當做新的特徵，這裡最高6次冪。程式碼如下

function out = mapFeature(X1, X2)
 % MAPFEATURE Feature mapping function to polynomial features
 %
 %   MAPFEATURE(X1, X2) maps the two input features
 %   to quadratic features used in the regularization exercise.
 %
 %   Returns a new feature array with more features, comprising of 
 %   X1, X2, X1.^2, X2.^2, X1*X2, X1*X2.^2, etc..
 %
 %   Inputs X1, X2 must be the same size
 %
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

？？？？？？？？？？看不懂

1. Cost function and gradient
   公式如下
   
   theta_0和之前的相同，其餘的引數在最後加了
   程式碼如下

hx = sigmoid(X * theta);
    J = (-1/m) * sum(y .* log(hx) + (1-y) .* log(1 - hx)) + ...
        (lambda/(2*m)) * (sum(theta .* theta) - theta(1)^2);

要注意減掉theta(1)的平方，因為公式中不包含 (theta_0)，所以需要去掉。

？？？？？？？？？？看不懂

梯度下降的演算法公式

計算程式碼

function [J, grad] = costFunctionReg(theta, X, y, lambda)
    m = length(y); % number of training examples
    grad = zeros(size(theta));

    hx = sigmoid(X * theta);
    J = (-1/m) * sum(y .* log(hx) + (1-y) .* log(1 - hx)) + ...
        (lambda/(2*m)) * (sum(theta .* theta) - theta(1)^2);

    grad(1) = (1/m) * sum((hx - y) .* X(:, 1));
    n = length(theta);
    for j = 2 : n
        grad(j) = (1/m) * sum((hx - y) .* X(:, j)) + (lambda / m) * theta(j);
    end
end

1. plotDecisionBoundary
   畫出非線性邊界的程式碼

else
    % Here is the grid range
    u = linspace(-1, 1.5, 50);
    v = linspace(-1, 1.5, 50);

    z = zeros(length(u), length(v));
    % Evaluate z = theta*x over the grid
    for i = 1:length(u)
        for j = 1:length(v)
            z(i,j) = mapFeature(u(i), v(j))*theta;
        end
    end
    z = z'; % important to transpose z before calling contour

    % Plot z = 0
    % Notice you need to specify the range [0, 0]
    contour(u, v, z, [0, 0], 'LineWidth', 2)
end

**

？？？？？？？？？？？？？看不懂

**
5. change regularization parameters

修改ex2_reg.m第90行可以修改lambda的大小，等於0對應不進行正則化。

作業的關鍵程式碼

plotData

pos = find(y==1);
neg = find(y==0);
plot(X(pos,1),X(pos,2),'k+','Linewidth',2,'MarkerSize',7);
plot(X(neg,1),X(neg,2),'ko','MarkerFaceColor','y','MarkerSize',7);

sigmoid

g=1./(1+exp(-z));

costFunction

z=X*theta;
h=sigmoid(z);
J=sum(-y.*log(h)-(1-y).*log(1-h))/m;

grad=(1/m)*(X.'*(h-y));

predict

h = sigmoid(X * theta); 
p = (h >= 0.5); 

costFunctionReg

h=sigmoid(X*theta);
J=sum(-y.*log(h)-(1-y).*log(1-h))/m + (lambda/(2*m))*sum(theta(2:end).^2);

grad(1)=(1/m)*(X(:,1).'*(h-y));
grad(2:end) =(1/m)*(X(:,2:end).'*(h-y))+lambda*theta(2:end)/m;