1: clear all; close all; clc

   2: x = load('ex4x.dat'); %每一行是一个样本

   3: y = load('ex4y.dat');

   4: [m, n] = size(x);

   5: sample_num = m;

   6: x = [ones(m, 1), x]; %x增加一维。因为前面文字或前几节说了

   7: % Plot the training data

   8: % Use different markers for positives and negatives

   9: figure;

  10: pos = find(y == 1); neg = find(y == 0);%pos 和neg 分别是 y元素=1和0的所在的位置序号组成的向量

  11: plot(x(pos, 2), x(pos,3), '+')%用+表示那些yi=1所对应的样本

  12: hold on

  13: plot(x(neg, 2), x(neg, 3), 'o')

  14: hold on

  15: xlabel('Exam 1 score')

  16: ylabel('Exam 2 score')

  17: itera_num=500;%迭代次数

  18: g = inline('1.0 ./ (1.0 + exp(-z))'); %这个就相当于制造了一个function g（z）=1.0 ./ (1.0 + exp(-z))

  19: plotstyle = {'b', 'r', 'g', 'k', 'b--', 'r--'};

  20: figure;%建立新的窗口

  21: alpha = [ 0.0009, 0.001,0.0011,0.0012,0.0013 ,0.0014 ];%下面就分别用这几个学习速率看看哪个更好

  22: for alpha_i = 1:length(alpha) %alpha_i是1,2，...6，表示的是学习速率向量和曲线格式向量的坐标：alpha(alpha_i)，plotstyle(alpha_i)

  23:     theta = zeros(n+1, 1);%thera表示样本Xi各个元素叠加的权重系数，这里以向量形式表示，且初始都为0，三维向量

  24:     J = zeros(itera_num, 1);%J是个100*1的向量，第n个元素代表第n次迭代cost function的值（下面用negtive 的对数似然函数，

  25:     %因为是negtive 的，所以是求得极小值）

  26:     for i = 1:itera_num %计算出某个学习速率alpha下迭代itera_num次数后的参数   

  27:         z = x * theta;%这个z是一个列向量，每一个元素是每一个样本Xi的线性叠加和，因为X是所有的样本，因此这里不是一个一个样本算的，

  28:         %而是所有样本一块算的，因此z是一个包含所有样本Xi的线性叠加和的向量。在公式中，是单个样本表示法，而在matlab中都是所有的样本一块来。

  29:         h = g(z);%这个h就是样本Xi所对应的yi=1时，映射的概率。如果一个样本Xi所对应的yi=0时，对应的映射概率写成1-h。

  30:         J(i) =(1/sample_num).*sum(-y.*log(h) - (1-y).*log(1-h));%损失函数的矢量表示法 这里Jtheta是个100*1的列向量。

  31:         grad = (1/sample_num).*x'*(h-y);%在程序中X是包括所有的X样本，公式中只能是每一个样本Xi，所以不能直接把公式照搬，所以要认真看看，代码中相应改变一下。

  33:         theta = theta - alpha(alpha_i).*grad;

  34:     end

35:     plot(0:itera_num-1, J(1:itera_num),char(plotstyle(alpha_i)),'LineWidth', 2)

%此处一定要通过char函数来转换因为包用（）索引后得到的还是包cell，

  36:     %所以才要用char函数转换，也可以用{}索引，这样就不用转换了。

  37:     %一个学习速率对应的图像画出来以后再画出下一个学习速率对应的图像。    

  38:     hold on

  39:     if(1 == alpha(alpha_i)) %通过实验发现alpha为0.0013 时效果最好，则此时的迭代后的theta值为所求的值

  40:         theta_best = theta;

  41:     end

  42: end

  43: legend('0.0009', '0.001','0.0011','0.0012','0.0013' ,'0.0014');%给每一个线段格式标注上

  44: xlabel('Number of iterations')

  45: ylabel('Cost function')

46: prob = g([1, 20, 80]*theta);

%把[1, 20, 80]*theta这个带入g（z）得到的就是在exam1=20、exam2=80的条件下，通过的概率（也就是Y=1）的概率是多少。

  47: %画出分界面

  48: % Only need 2 points to define a line, so choose two endpoints

  49: plot_x = [min(x(:,2))-2,  max(x(:,2))+2];%两个点就可以画出来直线，这里这么取x1坐标是让直接的视野更容易包含那些样本点。

50: plot_y = (-1./theta(3)).*(theta(2).*plot_x +theta(1));%分界面怎么画呢？问题也就是在x1，x2坐标图中找到

%那些将x1，x2带入1/(1+exp(-wx))后，使其值>0.5的（x1,x2）

  51: %坐标形成的区域，因为我们知道1/(1+exp(-wx))>0.5

52: %意味着该区域的（x1，x2）表示的成绩允许上大学的概率>0.5,那么其他的区域就是不被允许上大学，那么1/(1+exp(-wx))=0.5解出来的一个关

%于x1，x2的方程就是那个分界面。

  53: %我们解出来以后发现，这个方程是一个直线方程：w（2）x1+w（3）x2+w（1）=0

54: %注意我们不能因为这个分界面是直线，就认为logistic回归是一个线性分类器，注意logistic回归不是一个分类器，他没有分类的功能，

%这个logistic回归是用来

55: %预测概率的，这里具有分类功能是因为我们硬性规定了一个分类标准：把>0.5的归为一类，<0.5的归于另一类。这是一个很强的假设，

%因为本来我们可能预测了一个样本

56: %所属某个类别的概率是0.6，这是一个不怎么高的概率，但是我们还是把它预测为这个类别，只因为它>0.5.所以最后可能logistic回归加上这

%个假设以后形成的分类器

  57: %的分界面对样本分类效果不是很好，这不能怪logistic回归，因为logistic回归本质不是用来分类的，而是求的概率。

  58: figure;

  59: plot(x(pos, 2), x(pos,3), '+')%把分界面呈现在样本点上面，所以还要把样本点画出来。

  60: hold on

  61: plot(x(neg, 2), x(neg, 3), 'o')

  62: hold on

  63: plot(plot_x, plot_y)

  64: legend('Admitted', 'Not admitted', 'Decision Boundary')

  65: hold off

  66:  

y = load('ex4y.dat');[m, n] = size(x);% Add intercept term to x
x = [ones(m, 1), x]; % Plot the training data
% Use different markers for positives and negatives
figure
pos = find(y); neg = find(y == 0);%find是找到的一个向量，其结果是find函数括号值为真时的值的编号
plot(x(pos, 2), x(pos,3), '+')
hold on
plot(x(neg, 2), x(neg, 3), 'o')
hold on
xlabel('Exam 1 score')
ylabel('Exam 2 score')% Initialize fitting parameters
theta = zeros(n+1, 1);% Define the sigmoid function
g = inline('1.0 ./ (1.0 + exp(-z))'); % Newton's method
MAX_ITR = 7;
J = zeros(MAX_ITR, 1);for i = 1:MAX_ITR% Calculate the hypothesis functionz = x * theta;h = g(z);%是可以以向量形式给g（）函数输入的，不是因为所有函数都可以这么输，而是这个函数里面有个exp（）函数，%是可以的。z的每一维在exp（）中分别计算，对应h里的每一维。% Calculate gradient and hessian.% The formulas below are equivalent to the summation formulas% given in the lecture videos.grad = (1/m).*x' * (h-y);%H = (1/m).*x' * diag(h) * diag(1-h) * x;%因为在matlab中，是以所有样本一块作为X输入的，所以这里的h是一个%向量，里面的元素是每个样本映射的logistic回归模型的概率，这样两个对角矩阵乘出来，还是一个对角矩阵，每个对角%元是每个样本的hi*（1-hi）。这样做就是适应公式在matlab下的表现形式，是正确的，我验算了。% Calculate J (for testing convergence)J(i) =(1/m)*sum(-y.*log(h) - (1-y).*log(1-h));theta = theta - H\grad;%是这样的，matlab里 \ 是左除，/ 是右除。如果是数字计算，则左除和右除是等效的，例如 3/2 = 2\3。    %而对于矩阵运算，则二者不等效。矩阵除法在 matlab 里定义为矩阵求逆后相乘。例如 A的逆矩阵是 A1，则 B/A = B*A1，A\B = A1*B。    %矩阵乘法不满足交换律，因此需要有左右除法之分。矩阵求逆的命令是 inv矩阵乘法不满足交换律，因此需要有左右除法之分。

end
% Display theta
theta% Calculate the probability that a student with
% Score 20 on exam 1 and score 80 on exam 2
% will not be admitted
prob = 1 - g([1, 20, 80]*theta)%画出分界面
% Plot Newton's method result
% 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的计算公式见博客下面的评论。
plot_y = (-1./theta(3)).*(theta(2).*plot_x +theta(1));
plot(plot_x, plot_y)
legend('Admitted', 'Not admitted', 'Decision Boundary')
hold off% Plot J
figure
plot(0:MAX_ITR-1, J, 'o--', 'MarkerFaceColor', 'r', 'MarkerSize', 8)
xlabel('Iteration'); ylabel('J')
% Display J
J