分类的目标是预测离散值，获取输入向量x，并将其分给K个离散类Ck之一，其中k = 1，...，K。 对于离散的输入，我们从概率的角度考虑分类问题，如何从数据分布的简单假设中产生具有线性决策边界的模型。
主要讲了生成式模型和判别式模型。 举个例子：现在你要判断一群国际友人讲的是什么语言，有汉语、英语、西班牙语、法语等等。生成式模型就好比你去学习这些语言的模型，然后当你听到某个人讲什么语言的时候，你就可以判别出来。而判别式模型不一样，它是直接生成，相当于你找出来的是不同语言之间的的差别。听到某个人讲的语言，你能够通过判别来得到结果。生成式模型能推得判别式模型。反之，不能成立。
4.2概率生成模型
采用一种概率生成方法，通过训练数据得出类的条件密度p（x | Ck），以及类先验p（Ck），然后套用贝叶斯分类器使用这些数据计算后验概率p（Ck | x）。
考虑一个二分类问题：
类别1的后验概率可以表示为：
logistic sigmoid 函数又称挤压函数，因为它把整个实轴映射到一个有限的空间，是一个良好的阈值函数,连续，光滑,处处可导，
严格单调而且满足下面的对称性：
反函数是logit（分对数）：
考虑多分类的情况，我们有：
上面的式子称为归一化指数（ normalized exponential ）。其中：
归⼀化指数也被称为 softmax 函数，因为它表示“ max ”函数的一个平滑版本。 这是因为：如果对于所有的 j≠k，都有ak≫aj,
通常情况下softmax会被用在网络中的最后一层，用来进行最后的分类和归一化。
连续输入
我们假定类条件概率密度是高斯分布且所有的类别的协方差矩阵相同:
展开，类别 Ck 的类条件概率为:
D：变量维度； |Σ| ：是其行列式； Σ：协方差矩阵，描述各维变量之间的相关度。对于二维高斯分布，有：
首先考虑对于二分类高斯分布模型-生成模型：
其中：
（1）假设类概率的协方差矩阵相同，使得最后推导的结果中 x 的⼆次型消失了，决策边界是线性的。 （2）先验概率密度 p(Ck) 只出现在偏置中，说明先验概率密度变化的结果是决策边界平移，并不影响决策边界的形状。
上面的参数推导过程如下：
输入空间x是二维情况下的结果：
第一张图是两个不同的高斯分布，第二张图是对应的后验概率分布，由x的线性函数sigmoid函数给出。曲面颜色中，红色所占的比例为：p(C1|x),蓝色所占的比例：p(C2|x)=1-p(C1|x)。
再考虑K分类的高斯分布-生成模型：
这个和二分类问题很像，就是对应的等式是一个族，从k是1....K的范围。
ak(x)是x的线性函数，（各类别的协方差矩阵相同，使得⼆次项被消去）。最终的决策边界对应于最小错误分类率，会出现在后验概率最大的两个概率相等的位置。
不假设各个类别的协方差矩阵相同，允许每个类条件概率密度p(x | Ck)有自己的协方差矩阵Σk，那么之前⼆次项消去的现象不会出现，从而我们会得到x的2次函数。
第一张图是三个不同的高斯分布，绿色和红色的协方差矩阵是相同的，而蓝色的高斯分布是不同的。所以对应的后验概率分布图中，绿色和红色的决策边界是线性的，因为对应的二次项被消去了，而蓝色和其他的都是二次决策边界。
最大似然解
除了概率生成模型，最重要的还是最大似然估计的知识点。
极大似然估计：提供了一种给定观察数据来评估模型参数的方法，即：“模型已定，参数未知”。通过若干次试验，观察其结果，利用试验结果得到某个参数值能够使样本出现的概率为最大。将概率密度估计问题转化为参数估计问题。
前提：训练样本的分布能代表样本的真实分布。每个样本集的样本都是所谓独立同分布的随机变量。
方法：具体化条件概率密度p(x | Ck)参数化的函数，使用最大似然法确定参数的值以及类别先验概率p(Ck)。
两分类问题，每个类别都有一个高斯类条件概率密度，且协方差矩阵相同。我们假设我们有一个数据集 ：
其中n = 1, . . . , N。把先验概率记作p(C1) = π，p(C2) = 1 − π。 对于一个来自类别C1的数据点，我们有tn = 1，因此 对于类别C2的数据点，我们有tn = 0，因此：
对于类别C2的数据点，我们有tn = 0，因此：
对应的似然函数：
最大化似然函数的对数比较方便。首先考虑关于π的最大化。对数似然函数中与π相关的项为：
对变量π求导，令其为0：
这和我们所直观上判断的是一样的，应为相应的概率就是C1类的数量除以总的数量。
考虑关于µ1的最大化：
关于µ1、µ2的导数等于0，整理可得：
考虑协方差矩阵Σ的最大似然解：
求解结果：
使用高斯分布的最大似然解的标准结果，我们看到Σ = S，它表示对一个与两类都有关系的协方差矩阵求加权平均。
4.3 概率判别式模型
生成式模型的核心思想是通过最大似然法求出类概率密度的参数和类先验概率密度，然后用贝叶斯公式求出后验概率密度。是一种间接性的方法。 判别式模型将显示地使用⼀般的线性模型的函数形式，然后使用最大似然法直接确定它的参数。
固定基函数
基函数向量ϕ(x)对输入变量进行⼀个固定的非线性变换，最终的决策边界在特征空间ϕ中是线性的，而在原始x空间中是非线性决策边界， 基函数中的某⼀个通常设置为常数，例如ϕ0(x) = 1，使得对应的参数w0扮演偏置的作用。
给出了原始的输⼊空间(x1, x2)以及标记为红和蓝的数据点。这个空间中定义了两个“高斯”基函数ϕ1(x)和ϕ2(x)。右图给出了对应的特征空间(ϕ1, ϕ2)以及线性决策边界。决策边界由4.3.2节讨论的线性回归模型得到。对应的在原始空间中的非线性决策边界在左图中黑色曲线标记出。
logistic 回归
将二分类问题的后验概率写成作用在特征向量 ϕ 的线性函数上的 logistic sigmoid 函数的形式，即:
M 维特征空间 ϕ ： M 个可调节参数。
最大似然方法：2M+M(M+1)/2+1 个参数。2M:均值，协方差矩阵 对称，是一半的，1：先验概率
对于⼀个数据集（ϕn, tn），tn ∈ {0, 1}且ϕn = ϕ(xn)，并且n = 1, . . . , N：
最大似然方法无法区分某个解优于另⼀个解，并且在实际应用中哪个解被找到将会依赖于优化算法的选择和参数的初始化。即使数据点的数量多于参数，只要数据是线性可分的，这个问题就会出现。 解决办法：通过引入先验概率，然后寻找 ω 的 MAP 解。
迭代重加权最小平方求解：
得到了相应的梯度，我们要利用这个式子来进行权重的更新，误差函数是凸函数，因为Hassian矩阵是正定的，因此有一个唯一的最小值。此外，误差函数可以通过迭代重加权最小平方（IRLS）求出最小值，这种迭代方法基于Newton-Raphson迭代最优化框架：求解用到了迭代重加权最小平方：
对角矩阵是：
Hessian矩阵不再是常量，通过权矩阵R依赖于w。对应于误差函数不是二次函数的事实。 由0 < yn < 1，对于任意向量u都有 ，Hessian矩阵H是正定的。推出误差函数是w的凸函数，有唯一的最小值。
多类 logistic 回归：
在生成式模型中，后验概率由特征变量的线性函数的softmax变换给出：
我们考虑使用最大似然法直接确 定这个模型中的参数{wk}。使用“1-of-K”表达方式，属于类别Ck的特征向量ϕk的目标向量tn是一个二元向量，这个向量的第k个元素等于1，其余元素等于0。得到相应的似然函数：
先求出yk关于所有激活aj的导数。这些导数为：
求解结果：
梯度：误差(ynj − tnj )与基函数ϕn的乘积。这种梯度形式在线性模型的平方和误差函数以及logistic回归模型的误差函数中都出现过。 这个公式用于顺序算法。这种顺序算法中每次只出现一个模式，每个权向量都使用公式更新。
给出两种方法的对比：
最后，是相应的代码实现：
polynomial.py:
import itertools
import functools
import numpy as np
class PolynomialFeatures(object):
"""
polynomial features
transforms input array with polynomial features
Example
=======
x =
[[a, b],
[c, d]]
y = PolynomialFeatures(degree=2).transform(x)
y =
[[1, a, b, a^2, a * b, b^2],
[1, c, d, c^2, c * d, d^2]]
"""
def __init__(self, degree=2):
"""
construct polynomial features
Parameters
----------
degree : int
degree of polynomial
"""
assert isinstance(degree, int)
self.degree = degree
def transform(self, x):
"""
transforms input array with polynomial features
Parameters
----------
x : (sample_size, n) ndarray
input array
Returns
-------
output : (sample_size, 1 + nC1 + ... + nCd) ndarray
polynomial features
"""
if x.ndim == 1:
x = x[:, None]
x_t = x.transpose()
features = [np.ones(len(x))] # a list?
for degree in range(1, self.degree + 1):
for items in itertools.combinations_with_replacement(x_t, degree):
features.append(functools.reduce(lambda x, y: x * y, items)) # x*y
return np.asarray(features).transpose()classify.py:
import numpy as np
class Classifier(object):
"""
Base class for classifiers # frame
"""
def fit(self, X, t, **kwargs):
"""
estimate parameters given the training dataset
Parameters
----------
X : (sample_size, n_features) np.ndarray
training data input
t : (sample_size,) np.ndarray
training data target
"""
self._check_input(X)
self._check_target(t)
if hasattr(self, "_fit"):
self._fit(X, t, **kwargs)
else:
raise NotImplementedError
def classify(self, X, **kwargs):
"""
classify inputs
Parameters
----------
X : (sample_size, n_features) np.ndarray
samples to classify
Returns
-------
label : (sample_size,) np.ndarray
label index for each sample
"""
self._check_input(X)
if hasattr(self, "_classify"):
return self._classify(X, **kwargs)
else:
raise NotImplementedError
def proba(self, X, **kwargs):
"""
compute probability of input belonging each class
Parameters
----------
X : (sample_size, n_features) np.ndarray
samples to compute their probability
Returns
-------
proba : (sample_size, n_classes) np.ndarray
probability for each class
"""
self._check_input(X)
if hasattr(self, "_proba"):
return self._proba(X, **kwargs)
else:
raise NotImplementedError
def _check_input(self, X):
if not isinstance(X, np.ndarray):
raise ValueError("X(input) must be np.ndarray")
if X.ndim != 2:
raise ValueError("X(input) must be two dimensional array")
if hasattr(self, "n_features") and self.n_features != np.size(X, 1):
raise ValueError(
"mismatch in dimension 1 of X(input) (size {} is different from {})"
.format(np.size(X, 1), self.n_features)
)
def _check_target(self, t):
if not isinstance(t, np.ndarray):
raise ValueError("t(target) must be np.ndarray")
if t.ndim != 1:
raise ValueError("t(target) must be one dimensional array")
if t.dtype != np.int:
raise ValueError("dtype of t(target) must be np.int")
if (t < 0).any():
raise ValueError("t(target) must only has positive values")
def _check_binary(self, t):
if np.max(t) > 1 and np.min(t) >= 0:
raise ValueError("t(target) must only has 0 or 1")
def _check_binary_negative(self, t):
n_zeros = np.count_nonzero(t == 0)
n_ones = np.count_nonzero(t == 1)
if n_zeros + n_ones != t.size:
raise ValueError("t(target) must only has -1 or 1")
logistic_regression.py:
import numpy as np
from linear.classifier import Classifier # 引入Classifier框架
class LogisticRegressor(Classifier): # 参数传入 定义逻辑回归模型 离散的情况
"""
Logistic regression model
y = sigmoid(X @ w)
t ~ Bernoulli(t|y)
"""
def _fit(self, X, t, max_iter=100): # 输入样本 ， 0,1标签 ，最大迭代步数
self._check_binary(t)
w = np.zeros(np.size(X, 1)) # 初始化权重矩阵 X行
for _ in range(max_iter):
w_prev = np.copy(w) # 保存原先的权重信息 用来更新权重
y = self._sigmoid(X @ w) # sigmoid 特征向量@权重矩阵 输出y
grad = X.T @ (y - t) # 一阶导数
hessian = (X.T * y * (1 - y)) @ X # 二阶导数 Hessian矩阵
try:
w -= np.linalg.solve(hessian, grad) # 更新权重信息 学习率是1？
print(w)
except np.linalg.LinAlgError:
break
if np.allclose(w, w_prev): # 收敛到一定的精度 结束
break
self.w = w
def _sigmoid(self, a):
return np.tanh(a * 0.5) * 0.5 + 0.5 # tanh(a) = 2sigmoid(2a)-1
def _proba(self, X):
y = self._sigmoid(X @ self.w)
return y # 返回概率
def _classify(self, X, threshold=0.5):
proba = self._proba(X)
label = (proba > threshold).astype(np.int)
return label # 返回分类结果 大于0.5 1 否则 0
class BayesianLogisticRegressor(LogisticRegressor): # 连续的情况 分布服从高斯分布
"""
Logistic regression model
w ~ Gaussian(0, a^(-1)I)
y = sigmoid(X @ w)
t ~ Bernoulli(t|y)
"""
def __init__(self, alpha=1.): # alpha = 1
self.alpha = alpha
def _fit(self, X, t, max_iter=100): # 输入特征向量 样本标签 最大迭代步数
self._check_binary(t) # 是否是 0 1 情况
w = np.zeros(np.size(X, 1))
eye = np.eye(np.size(X, 1))
self.w_mean = np.copy(w)
self.w_precision = self.alpha * eye
for _ in range(max_iter):
w_prev = np.copy(w)
y = self._sigmoid(X @ w)
grad = X.T @ (y - t) + self.w_precision @ (w - self.w_mean)
hessian = (X.T * y * (1 - y)) @ X + self.w_precision
try:
w -= np.linalg.solve(hessian, grad)
except np.linalg.LinAlgError:
break
if np.allclose(w, w_prev):
break
self.w_mean = w
self.w_precision = hessian
def _proba(self, X):
mu_a = X @ self.w_mean
var_a = np.sum(np.linalg.solve(self.w_precision, X.T).T * X, axis=1)
y = self._sigmoid(mu_a / np.sqrt(1 + np.pi * var_a / 8))
return y
softmax_regressor.py:
import numpy as np
from linear.classifier import Classifier
class SoftmaxRegressor(Classifier):
"""
Softmax regression model
aka multinomial logistic regression,
multiclass logistic regression, or maximum entropy classifier.
y = softmax(X @ W)
t ~ Categorical(t|y)
"""
def _fit(self, X, t, max_iter=100, learning_rate=0.1):
self.n_classes = np.max(t) + 1
T = np.eye(self.n_classes)[t]
W = np.zeros((np.size(X, 1), self.n_classes))
for _ in range(max_iter):
W_prev = np.copy(W)
y = self._softmax(X @ W)
grad = X.T @ (y - T)
W -= learning_rate * grad
if np.allclose(W, W_prev):
break
self.W = W
def _softmax(self, a):
a_max = np.max(a, axis=-1, keepdims=True)
exp_a = np.exp(a - a_max)
return exp_a / np.sum(exp_a, axis=-1, keepdims=True)
def _proba(self, X):
y = self._softmax(X @ self.W)
return y
def _classify(self, X):
proba = self._proba(X)
label = np.argmax(proba, axis=-1)
return label
测试代码test.py：
# -*- coding:UTF-8 -*-
import numpy as np
import matplotlib.pyplot as plt
from features import PolynomialFeatures # 多项式值 x-->1，x1,x2,x1*x2,...
from linear import LogisticRegressor, LeastSquaresClassifier # 用到了这个
np.random.seed(1234)
def create_data(add_outliers=False, add_class=False):
x0 = np.random.normal(size=50).reshape(-1, 2) - 1 # 正态分布 50 未知行数 2列 -1
x1 = np.random.normal(size=50).reshape(-1, 2) + 1 # 正态分布 50 未知行数 2列 +1
if add_outliers:
x_1 = np.random.normal(size=10).reshape(-1, 2) + np.array([5., 10.])
return np.concatenate([x0, x1, x_1]), np.concatenate([np.zeros(25), np.ones(30)]).astype(np.int)
if add_class:
x2 = np.random.normal(size=50).reshape(-1, 2) + 3
return np.concatenate([x0, x1, x2]), np.concatenate(
[np.zeros(25), np.ones(25), 2 + np.zeros(25)].astype(np.int))
return np.concatenate([x0, x1]), np.concatenate([np.zeros(25), np.ones(25)]).astype(np.int) # 0 1
x_train, y_train = create_data(add_outliers=True) # x_train(x0,x1) y_train(00...000,11...111)
# x_train, y_train = create_data(add_outliers=False) # x_train(x0,x1) y_train(00...000,11...111)
x1_test, x2_test = np.meshgrid(np.linspace(-5, 15, 100), np.linspace(-5, 15, 100)) # 产生测试集数据,-5到5 100个
x_test = np.array([x1_test, x2_test]).reshape(2, -1).T # 测试集 reshape T
feature = PolynomialFeatures(1) # 线性多项式
X_train = feature.transform(x_train) #
# print(X_train)
X_test = feature.transform(x_test)
# print(X_test)
# -----------------logistic regression-----------------------#
model_logisticregression = LogisticRegressor()
model_logisticregression.fit(X_train, y_train)
y_lr = model_logisticregression.proba(X_test)
plt.scatter(x_train[:, 0], x_train[:, 1], c=y_train) # 绘制数据散点图 x,y 图像是0 1两种颜色
plt.contourf(x1_test, x2_test, y_lr.reshape(100, 100), alpha=0.2, levels=np.linspace(0, 1, 5)) # 模型输出的概率
plt.colorbar() #
plt.ylim(-5, 15) # 范围
plt.xlim(-5, 15)
plt.title("LogisticRegressor")
plt.gca().set_aspect('equal', adjustable='box')
plt.show()
# ------------------LeastSquaresClassifier---------------------#
model_least_squares = LeastSquaresClassifier()
model_least_squares.fit(X_train, y_train)
y_ls = model_least_squares.classify(X_test)
plt.scatter(x_train[:, 0], x_train[:, 1], c=y_train)
plt.contourf(x1_test, x2_test, y_ls.reshape(100, 100), alpha=0.2, levels=np.linspace(0, 1, 5))
plt.xlim(-5, 15)
plt.ylim(-5, 15)
plt.colorbar() #
plt.gca().set_aspect('equal', adjustable='box')
plt.title("LeastSquaresClassifier")
plt.show()
结果：