ML Note 1.2 - Classification

对于离散型随机变量的学习问题称为 classification.

Generative Learning Algorithms

生成式模型中，我们首先对 class-conditioned probability density $p(x|y)$ 和 prior $p(y)$ 进行建模。这一步可以对联合分布应用 MLE 来完成
$$l(\theta )=\log \prod _{i}p(x^{(i)},y^{(i)})$$

而后我们使用贝叶斯公式计算后验概率
$$p(y|x)=\frac{p(x|y)p(y)}{p(x)}\propto p(x|y)p(y)$$

一般来说，预测时只需要最小化错误率即可
$$h(x)=\underset{y}{\max }p(x|y)p(y)$$

有时一些不正确的预测结果可能会导致较大损失。为了尽可能减少预测错误带来的损失，可以构造决策表 $\lambda (h(x),y)$ 表示当实际类别为 $y$ 而预测类别为 $h(x)$ 时的损失。定义条件期望损失
$$R(\alpha |x)=\sum _i\lambda (\alpha ,i)p(x|i)p(i)$$

则预测可以表示成
$$h(x)=\underset{\alpha }{\min }R(\alpha |x)$$

Gaussian Discriminant Analysis

假设对于二分类问题
$$\begin{array}{rcl}y& \sim & Bern(\varphi )\\ x|y=1& \sim & MVN({\vec{\mu }}_1,\Sigma )\\ x|y=0& \sim & MVN({\vec{\mu }}_0,\Sigma )\end{array}$$

MLE 给出
$$\begin{array}{rcl}\varphi & =& \frac{1}{m}\sum _{i=1}^m{y}^{(i)}\\ {\vec{\mu }}_0& =& \frac{\sum _{i=1}^m\left(1-y^{(i)}\right)x^{(i)}}{\sum _{i=1}^m\left(1-y^{(i)}\right)}\\ {\vec{\mu }}_1& =& \frac{\sum _{i=1}^m{y}^{(i)}{x}^{(i)}}{\sum _{i=1}^m{y}^{(i)}}\\ \Sigma & =& \frac{1}{m}\sum _{i=1}^m\left(x^{(i)}-{\mu }_{y^{(i)}}\right){\left(x^{(i)}-{\mu }_{y^{(i)}}\right)}^T\end{array}$$

GDA 的后验概率可以写成如下形式
$$p(y=1|x;\varphi ,\Sigma ,{\mu }_0,{\mu }_1)=\frac{1}{1+\exp (-{\theta }_{\varphi ,\Sigma ,{\mu }_0,{\mu }_1}^{T}x)}$$

尽管形式上与逻辑回归相同，但一般来说二者的决策边界并不相同。事实上，任何指数族分布 $x|y\sim ExpFamily(\eta )$ 都将推出对数后验。因此逻辑回归有更强的 robust 性能，而且对模型的假设较少。而 GDA 假设数据来自正态分布，因此收敛速度较快。

Naive Bayes Classifier

如果特征的各个分量 $x_i$ 之间条件独立，即满足 NB assumption
$$p(x_1,\dots ,x_n|y)=\prod _{i=1}^{n}p(x_i|y)$$

则可以使用朴素贝叶斯分类器进行建模。考虑如下的问题

判别一封邮件是否是垃圾邮件。假设我们有一本字典 $V$ 可以将单词映射为整数 { 1 , … , n } \{1, \dots, n\} {1,…,n}.

取特征 x ∈ { 0 , 1 } n x\in \{0,1\}^n x∈{0,1}n 表示字典中第 $i$ 个单词是否出现。应用 MLE 可以得到参数
$$\begin{array}{rcccl}{\varphi }_{i|y=1}& \equiv & p(x_i=1|y=1)& =& \frac{\sum _{j=1}^m{x}_i^{(j)}{y}^{(j)}}{\sum _{j=1}^m{y}^{(j)}}\\ {\varphi }_{i|y=0}& \equiv & p(x_i=1|y=0)& =& \frac{\sum _{j=1}^m{x}_i^{(j)}(1-y^{(j)})}{\sum _{j=1}^m(1-y^{(j)})}\\ {\varphi }_y& \equiv & p(y=1)& =& \frac{\sum _{j=1}^m{y}^{(j)}}{m}\end{array}$$

上述模型被称为 multi-variate Bernoulli event model 因为特征 $x$ 是一个多元伯努利变量。另一种可能的模型被称为 multinomial event model 其中特征 x ∈ { 1 , … , n } l x \in \{1,\dots,n\}^l x∈{1,…,n}l 是对于邮件的一个翻译。其中 $l$ 表示邮件的单词数。应用 MLE 可以得到参数
ϕ i ∣ y = 1 ≡ p ( x j = i ∣ y = 1 ) = ∑ i = 1 m ∑ j = 1 l i y ( i ) I { x j ( i ) = k } ∑ i = 1 m y ( i ) l i ϕ i ∣ y = 0 ≡ p ( x j = i ∣ y = 0 ) = ∑ i = 1 m ∑ j = 1 l i ( 1 − y ( i ) ) I { x j ( i ) = k } ∑ i = 1 m ( 1 − y ( i ) ) l i ϕ y ≡ p ( y = 1 ) = ∑ i = 1 m y ( i ) m \begin{array}{rcccl} \phi_{i|y=1} &\equiv& p(x_j=i|y=1) &=& \frac{\sum\limits_{i=1}^m \sum\limits_{j=1}^{l_i}y^{(i)}I\{x_j^{(i)}=k\}}{\sum\limits_{i=1}^my^{(i)}l_i}\\ \phi_{i|y=0} &\equiv& p(x_j=i|y=0) &=& \frac{\sum\limits_{i=1}^m\sum\limits_{j=1}^{l_i}(1-y^{(i)})I\{x_j^{(i)}=k\}}{\sum\limits_{i=1}^m(1-y^{(i)})l_i}\\ \phi_y &\equiv& p(y=1) &=& \frac{\sum\limits_{i=1}^my^{(i)}}{m} \end{array} ϕi∣y=1​ϕi∣y=0​ϕy​​≡≡≡​p(xj​=i∣y=1)p(xj​=i∣y=0)p(y=1)​===​i=1∑m​y(i)li​i=1∑m​j=1∑li​​y(i)I{xj(i)​=k}​i=1∑m​(1−y(i))li​i=1∑m​j=1∑li​​(1−y(i))I{xj(i)​=k}​mi=1∑m​y(i)​​

上述算法的问题在于，如果一个单词在所有垃圾邮件样本中都没有出现，则待测样本中只要出现这个单词就一定不会被判为垃圾邮件。这个问题可以被 Laplace smoothing 解决，其思想是给每个单词一个默认的出现概率。设随机变量 $z\sim Multinomial({\varphi }_1,{\varphi }_2,\dots ,{\varphi }_{k-1})$ 的 MLE 结果为 ${\varphi }_j=p/q$ 。因为 $z$ 从 $k$ 个候选中取值，因此每个值的默认出现概率应为 $1/k$ 。修改后的参数估计值为
$${\varphi }_j=\frac{p+1}{q+k}$$

应该注意到，上面的算法只对离散型随机变量 $x$ 有效。如果希望对连续型随机变量应用朴素贝叶斯，可以将其 discretize 成为离散型随机变量。

Discriminative

判别式模型中，我们直接对后验概率 $p(y|x)$ 建模，并以 $p(y|x)=0.5$ 作为决策边界。

Softmax Regression

对于 $k$ 类决策问题，假设目标变量
$$y|x;\theta \sim Multinomial(1,\vec{\varphi })$$

其中每一类的出现概率
$$\vec{\varphi }=({\varphi }_1,{\varphi }_2,\dots ,{\varphi }_{k-1})$$

一般不将 ${\varphi }_k$ 作为参数之一，这是因为归一性要求
$${\varphi }_k=1-\sum _{i=1}^{k-1}{\varphi }_i$$

现在我们希望证明多项分布属于指数族。为了使
$$h(x)=E[T(y)]=\left[\begin{array}{c}{\varphi }_1\\ {\varphi }_2\\ ⋮\\ {\varphi }_{k-1}\end{array}\right]$$

构造
T ( y ) = [ I { y = 1 } I { y = 2 } ⋮ I { y = k − 1 } ] T(y) = \left[\begin{array}{c} I\{y = 1\}\\ I\{y = 2\}\\ \vdots\\ I\{y = k-1\}\\ \end{array}\right] T(y)=⎣⎢⎢⎢⎡​I{y=1}I{y=2}⋮I{y=k−1}​⎦⎥⎥⎥⎤​

则可以观察到
I { y = k } = 1 − ∑ i = 1 k − 1 T ( y ) I\{y = k\} = 1 - \sum\limits_{i = 1}^{k - 1} T(y) I{y=k}=1−i=1∑k−1​T(y)

因为多项分布满足
p ( y ; ϕ ) = ∏ i = 1 k ϕ i I { y = i } = exp ⁡ ( ∑ i = 1 k I { y = i } log ⁡ ϕ i ) = exp ⁡ ( ∑ i = 1 k − 1 T i ( y ) log ⁡ ϕ i + ( 1 − ∑ i = 1 k − 1 T i ( y ) ) log ⁡ ϕ k ) = exp ⁡ ( ∑ i = 1 k − 1 T i ( y ) log ⁡ ϕ i ϕ k + log ⁡ ϕ k ) \begin{array}{rcl} p(y;\phi) &=& \prod\limits_{i=1}^k\phi_i^{I\{y = i\}}\\ &=& \exp\left(\sum\limits_{i = 1}^k I\{y = i\}\log\phi_i\right)\\ &=& \exp\left(\sum\limits_{i = 1}^{k - 1} T_i(y)\log\phi_i + \left(1 - \sum\limits_{i = 1}^{k - 1} T_i(y)\right)\log\phi_k\right)\\ &=& \exp\left(\sum\limits_{i = 1}^{k - 1} T_i(y)\log\frac{\phi_i}{\phi_k} + \log\phi_k\right) \end{array} p(y;ϕ)​====​i=1∏k​ϕiI{y=i}​exp(i=1∑k​I{y=i}logϕi​)exp(i=1∑k−1​Ti​(y)logϕi​+(1−i=1∑k−1​Ti​(y))logϕk​)exp(i=1∑k−1​Ti​(y)logϕk​ϕi​​+logϕk​)​

所以指数族分布参数为
$$\begin{array}{rcl}\eta & =& \left[\begin{array}{c}\log ({\varphi }_1/{\varphi }_k)\\ \log ({\varphi }_2/{\varphi }_k)\\ ⋮\\ \log ({\varphi }_{k-1}/{\varphi }_k)\end{array}\right]\\ a(\eta )& =& -log({\varphi }_k)\\ b(y)& =& 1\end{array}$$

定义 ${\eta }_k\equiv \log \frac{{\varphi }_k}{{\varphi }_k}=0$ 则可以证明
$$h(x)=\left[\begin{array}{c}\frac{\exp {\eta }_1}{\sum _{j=1}^k\exp {\eta }_j}\\ \frac{\exp {\eta }_2}{\sum _{j=1}^k\exp {\eta }_j}\\ ⋮\\ \frac{\exp {\eta }_{k-1}}{\sum _{j=1}^k\exp {\eta }_j}\end{array}\right]=\left[\begin{array}{c}\frac{\exp {\theta }_1^{T}x}{\sum _{j=1}^k\exp {\theta }_j^{T}x}\\ \frac{\exp {\theta }_2^{T}x}{\sum _{j=1}^k\exp {\theta }_j^{T}x}\\ ⋮\\ \frac{\exp {\theta }_{k-1}^{T}x}{\sum _{j=1}^k\exp {\theta }_j^{T}x}\end{array}\right]$$

Logistic Regression

对于二分类问题，假设目标变量
$$y|x;\theta \sim Bern(p)$$

为 softmax regression 在 $k=2$ 时的特殊情况。因此可以得出
$$h(x)=\frac{1}{1+e^{-{\theta }^{T}x}}$$

这样的函数被称为 logistic 或 sigmoid
$$g(z)=\frac{1}{1+e^{-z}}$$

所以该模型称为 logistic regression. 应用 MLE 得到
$$\begin{array}{rcl}L(\theta )& =& \prod _{i=1}^{m}g{\left({\theta }^T{x}^{(i)}\right)}^{y^{(i)}}{\left(1-g\left({\theta }^T{x}^{(i)}\right)\right)}^{1-y^{(i)}}\\ l(\theta )& =& \sum _{i=1}^m{y}^{(i)}\ln g\left({\theta }^T{x}^{(i)}\right)+\left(1-y^{(i)}\right)\ln \left(1-g\left({\theta }^T{x}^{(i)}\right)\right)\end{array}$$

根据 sigmoid 函数的对称性
$$g(-z)=1-g(z)$$

将上式化为
$$\begin{array}{rcl}l(\theta )& =& \sum _{i=1}^m{y}^{(i)}\ln g\left({\theta }^T{x}^{(i)}\right)+\left(1-y^{(i)}\right)\ln g\left(-{\theta }^T{x}^{(i)}\right)\\ & =& \sum _{i=1}^m{y}^{(i)}{\theta }^T{x}^{(i)}+\ln g\left(-{\theta }^T{x}^{(i)}\right)\end{array}$$

因为 sigmoid 函数的导数满足
$$g^{\prime }(z)=g(z)(1-g(z))$$

所以
$$\begin{array}{rcl}{\nabla }_{\theta }l& =& \sum _{i=1}^m{y}^{(i)}{x}^{(i)}+\frac{1}{g\left(-{\theta }^T{x}^{(i)}\right)}\cdot g\left(-{\theta }^T{x}^{(i)}\right)\left(1-g\left(-{\theta }^T{x}^{(i)}\right)\right)\left(-x^{(i)}\right)\\ & =& \sum _{i=1}^m{y}^{(i)}{x}^{(i)}-g\left({\theta }^T{x}^{(i)}\right)x^{(i)}\\ & =& \sum _{i=1}^m\left(y^{(i)}-g\left({\theta }^T{x}^{(i)}\right)\right)x^{(i)}\end{array}$$

由此就可以通过梯度下降来求解了。

---

Appendix

Logistic Regression

main.m

%% Initialization
clear ; close all; clc

%% Load Data
%  The first two columns contains the X values and the third column
%  contains the label (y).

data = load('ex2data2.txt');
X = data(:, [1, 2]); y = data(:, 3);

% Note that mapFeature also adds a column of ones for us, so the intercept
% term is handled
X = mapFeature(X(:,1), X(:,2));

%% Regularization and Accuracies

initial_theta = zeros(size(X, 2), 1);
lambda = 1;

% Set Options
options = optimset('GradObj', 'on', 'MaxIter', 400);

% Optimize
[theta, J, exit_flag] = ...
	fminunc(@(t)(costFunctionReg(t, X, y, lambda)), initial_theta, options);

%% Plot Boundary

% Plot Data
figure; hold on;
pos = (y == 1);
neg = (y == 0);

plot(X(pos, 2), X(pos, 3), 'k+');
plot(X(neg, 2), X(neg, 3), 'ko');
hold on

% 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)

title(sprintf('lambda = %g', lambda))

% Labels and Legend
xlabel('Microchip Test 1')
ylabel('Microchip Test 2')

legend('y = 1', 'y = 0', 'Decision boundary')
hold off;

%% Compute accuracy on our training set
p = double(logsig(X * theta) >= 0.5);

fprintf('Train Accuracy: %f\n', mean(double(p == y)) * 100);
fprintf('Expected accuracy (with lambda = 1): 83.1 (approx)\n');

costFunctionReg.m

function [J, grad] = costFunctionReg(theta, X, y, lambda)
%COSTFUNCTIONREG Compute cost and gradient for logistic regression with regularization
%   J = COSTFUNCTIONREG(theta, X, y, lambda) computes the cost of using
%   theta as the parameter for regularized logistic regression and the
%   gradient of the cost w.r.t. to the parameters. 

% Initialize some useful values
m = length(y); % number of training examples

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the cost of a particular choice of theta.
%               You should set J to the cost.
%               Compute the partial derivatives and set grad to the partial
%               derivatives of the cost w.r.t. each parameter in theta

J = -1 / m * (sum(y .* log(logsig(X * theta)) + (1 - y) .* log(1 - logsig(X * theta))) - lambda / 2 * norm(theta(2:end))^2);
grad = 1 / m * (X' * (logsig(X * theta) - y) + lambda * [0; theta(2:end)]);

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

end

mapFeature.m

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