现在的位置: 首页 > 烟台股票配资公司 > 正文
A- A+
广义线性模型与逻辑回归-四川大熊猫栖息地

零、概述

本文内容包括:线性模型、广义线性模型、逻辑回归、牛顿法以及 IRLS 。前导知识包括:矩阵、向量的计算,多元函数的梯度、赫森矩阵以及二阶泰勒展开,概率论。


一、线性模型

假设样本由 维特征向量 和标量响应 组成。 和 服从联合分布 。把确定 的前提下 的条件分布记作 。如果假设对于 , 的随机噪声是加性的,如下式:

则当 时 本质上由 确定,其所有不确定性来源于 。如果进一步假设 满足高斯分布:

那么这就是加性高斯误差。在这种设定下分布 的条件期望,也是概率密度最大的值就是:

把该条件期望当作 时对 的预测值。至于 则是问题固有的、无法通过模型优化去除的误差。所谓线性模型就是用线性函数来建模 :

最后一个等号后面的 是 维向量 , 是 维向量 。 前面加一维常量 1 是为了把截距 纳入向量表达式中。训练集( )是 个特征向量及其标量响应: 。如果问题是一个回归问题,且选用平方误差作为损失函数:

令 在训练集上最小化平均损失:

最后一个等号后面的 是所有 组成的 维向量。 是 矩阵,其第 行是 。损失函数作为 的函数是凸的,存在全局唯一最小点,也就是唯一的驻点。通过求导并置导向量为 0 再解方程可得到:

这就是要寻找的最优参数——线性回归的最小二乘解。当 时对 个训练样本的预测值是:

是 的列的线性组合,系数向量是 。也就是说 属于 的列空间。并且由于 最优化了平均平方误差,所以 是 的列空间中与向量 欧氏距离最小的向量,即 是 在 列空间上的投影。同时 也是 的一个线性变换,变换矩阵是 。这也是该模型被称为线性模型的一个原因。 的对角线元素的和,即 的迹为:

上式第二个等号利用了 的迹与 的迹相等。把矩阵乘积的迹的表达式按元素展开很容易证明这一点。 的迹等于自由参数的个数。所以把 的迹定义为线性模型的自由度是合适的。如果为模型加上正则化,那么 就受到了一种约束,它在参数空间中的变化就没有那么“自由”了。正则化的线性模型的自由度就不应该是单纯的参数个数。其变换矩阵的迹小于 ,用迹度量正则化线性模型的自由度仍是恰当的(参见本专栏文章 理解机器学习中的 L2 正则化)。


二、广义线性模型和逻辑回归

把 在 时的条件期望简记作:

如果引入一个函数 ,令:

称作链接函数 。如果 是恒等函数 ,则 [2.2] 就退化为 [1.4] 。所以线性模型是 [2.2] 在链接函数为恒等函数时的特例。[2.2] 就是广义线性模型。链接函数 采用不同的形式就产生不同的模型。

考虑二分类问题。样本输出 是类别的标记:0 表示负例;1 表示正例。把 时样本属于正例的概率记为 ,则有:

最后的 是 的简写。令链接函数 是 函数:

函数图形为:

图 2.1 logit 函数图形

把 函数和 [2.3] 代入 [2.2] 中有:

也就是将正例与负例概率之比的对数(称作 )建模为线性。从这里也可看出,不论最终将判断正例的概率阈值定在多少,模型判断正例与负例的分界线是 ,为一条直线。现在将 [2.5] 变形得到:

这里的记法 强调一下 受 影响。[2.6] 就是逻辑回归模型。可见逻辑回归属于广义线性模型。分类问题与回归问题表面上看是不同的问题,但在深层次上它们是统一的——都是用损失函数衡量拟合的损失,调整参数尽量降低这个损失。

对于回归问题,可以用平方误差 [1.5] 衡量拟合的损失。当模型输出不等于真实值时应该加以惩罚。平方误差在模型输出不等于真实值时大于 0 ,且相差越大则平方误差越大。但是平方误差随着模型预测与真实值差距的增大呈二次增长,对于大差距惩罚得太过。如果存在离群样本,比如测量错误或标注错误,平方误差易受干扰,并不健壮。存在其他惩罚函数以缓和这个问题。分类问题的损失函数又应该有另外的形式。

对于分类问题,现在令 取 1 或 -1 标记正负例。模型 是连续函数。如果当 时将样本判为正例;当 时判为负例。那么当 和 同号,即 时模型分类正确;当 时模型分类错误。注意当 只能取 1 或 -1 时,平方误差可以视作 的函数:

其图形为(蓝色曲线):

图 2.2 三种损失函数的图形。黑色虚线是错分类惩罚:分错惩罚为 1.0;分对惩罚为 0.0 。

观察蓝色曲线。当 时分类错误。如果把 的绝对值作为分类置信程度的话,负的 越小错误越严重。蓝色线条在这种情况下升高,惩罚这种现象。当 时分类是正确的。且 越大则分得越对。这时惩罚应该减小。但是平方误差的行为是先减小,到大于 1.0 后反而增大。平方误差是评价回归质量的。回归问题小也不行大也不行。但是分类问题性质不同,平方误差不再恰当。

对于二分类问题,一个可采用的损失函数是 。如果 取 1 和 0 标记正负例。令 时样本为正例的概率是 ,该损失函数为:

它是“对数似然”的相反数。将 [2.8] 取相反就是对数似然。可以到看当 时对数似然为 ;当 时对数似然为 。 也是交叉熵 。交叉熵衡量两个分布的相似程度。其中一个分布是模型的预测分布。另一个分布是观察到的该样本类别的分布。对数似然越大越好。最小化 就是最大化对数似然,最小化交叉熵。

如果用 来标记正负例,则 在正负例时取值分别是 1 和 0 。将 和式 [2.6] 带入 [2.8] ,令 ,有:

第二步推导有些跳。记住 的值只有两种可能 1 或 -1 。推导不下去时把式中的 替换成 1 和 -1 分别继续推导。最后两条路合成一个最终的结果就是 [2.9] 。 的图形是图 2.2 中的橙色曲线:当 时增大,但没有误差平方增大得那么快。所以它对于离群样本更健壮。当 时减小,且随着 增大越来越小。这说明对于分类问题来说 的惩罚行为是恰当的。

以下继续用 标记正负例。采用 作为损失函数。逻辑归回就是寻找 在训练集上最小化平均损失:

当找到了 ,一个新的样本 属于正例的概率就是:

可以根据 是否大于 0.5 (或者其他阈值)判断 是否属于正例。从这个式子也可以看出逻辑回归模型相当于只有一个神经元的神经网络,该神经元激活函数是 。 不是凸函数,其最小值没有解析解,需要通过最优化算法来寻找 ,比如梯度下降法(参见本专栏文章 神经网络之梯度下降与反向传播(上))。本文介绍基于目标函数二阶性质的牛顿法来寻找逻辑回归的解。


三、牛顿法与迭代重加权最小二乘(IRLS)

如果一个函数 是二阶可导的,在任一个数据点 附近可以将它二阶泰勒展开:

是 函数。 和 都是 维向量。 是 在 处的梯度,是一个 维向量。 是 在 处的二阶导矩阵——赫森矩阵,是一个 的矩阵。牛顿法的原理就是在一个点将函数二阶泰勒展开,展开后得到一个凸二次函数,是原函数在该点附近的二阶近似。将该二阶近似的全局最小点作为迭代的下一个点。如图所示:

图 3.1 牛顿法一步迭代示意图

二次函数在全局最小点梯度为 0 。将 [3.1] 约等号右侧求导后置 0 并解方程求得全局最小点是:

[3.2] 就是牛顿法的迭代公式。如果 本身就是一个二次函数(例如平方误差 [1.6] ) ,则它在任一点的二阶泰勒展开就是它本身。这种情况下牛顿法只用一步迭代就可以到达全局最小点。最小二乘解 [1.7] 就是牛顿法一步达到的 。

现在应用牛顿法来寻找逻辑回归模型的最优参数。为了方便去掉了 系数,这不影响寻找最小值:

将 [3.3] 式对 求导(并转置)得到 在 的梯度:

将梯度求导得到 在 的赫森矩阵:

[3.5] 有些跳步,因为中间过程的式子不好编辑。推导时注意 对 的导数是 。[3.4] 和 [3.5] 可以整理成矩阵形式:

为 矩阵,第 行为 。 和 分别是由 和 组成 维向量。 是 对角矩阵,对角线元素是 。将 [3.6] 代入牛顿法迭代公式 [3.2] 得到 :

[3.7] 就是用牛顿法训练逻辑回归模型的迭代公式。训练过程直到 收敛时停止。继续将 [3.7] 变形:

[3.8] 中的 。可以看出 是以 为特征以 为响应的加权最小二乘解。 被称为工作响应。 每一次迭代都有变化,是迭代过程中的临时响应。权值矩阵 也随着每次迭代而变化。故用 [3.8] 求逻辑回归最优参数又被称为迭代重加权最小二乘,Iteratively Reweighted Least Squares,IRLS 。


四、实现

以下用 python 实现逻辑回归。实现了梯度下降、 IRLS 和通用牛顿法三种优化算法。支持 mini batch 。通用牛顿法是指采用 [3.7] 式对 进行更新。牛顿法和梯度下降法可以通过 控制学习速率并支持学习率衰减。梯度下降支持 L2 正则化和冲量。(请见专栏文章)。

import numpy as np


class LogisticRegression:
    def __init__(self, method="Gradient", eta=0.1, threshold=1e-5, max_epochs=10, regularization=0.1, minibatch_size=5,
                 momentum=0.9, decay_power=0.5,
                 verbose=False):

        self.method = str(method)  # IRLS, Newton or Gradient
        self.eta = float(eta)  # learning rate for Newton and Gradient method
        self.threshold = float(threshold)  # stopping loss threshold
        self.max_epochs = int(max_epochs)  # max number of iterations
        self.regularization = float(regularization)  # L2 regularization strength, only for Gradient Descent
        self.minibatch_size = int(minibatch_size)  # minibatch size
        self.decay_power = float(decay_power)  # learning rate decaying power for gradient descent and newton
        self.momentum = float(momentum)  # gradient descent momentum
        self.verbose = verbose

        self.beta = None

    def fit(self, x, y):
        x = np.mat(np.c_[[1.0] * x.shape[0], x])
        y = np.mat(y).T
        I = np.eye(x.shape[1])
        self.beta = np.mat(np.random.random((x.shape[1], 1)) / 100)
        accu_gradient = np.mat(np.zeros(self.beta.shape))

        epochs = 0
        iterations = 0
        start = 0
        train_set_size = x.shape[0]
        loss = []

        effective_eta = self.eta
        while True:

            end = start + self.minibatch_size
            minibatch_x = x[start:end]
            minibatch_y = y[start:end]
            start = (start + self.minibatch_size) % train_set_size

            linear_combination = -np.matmul(minibatch_x, self.beta)
            linear_combination[linear_combination > 1e2] = 1e2  # prevent overflow
            p = 1.0 / (1.0 + np.power(np.e, linear_combination))
            p[p >= 1.0] = 1.0 - 1e-10
            p[p <= 0.0] = 1e-10  # truncate the probabilities to prevent W from being singular
            w = np.mat(np.diag(np.multiply(p, 1.0 - p).A1))

            if self.method == "IRLS":
                z = minibatch_x * self.beta + w.I * (minibatch_y - p)
                self.beta = (minibatch_x.T * w * minibatch_x + 1e-15 * I).I * minibatch_x.T * w * z
            elif self.method == "Newton":
                self.beta = self.beta + effective_eta * (minibatch_x.T * w * minibatch_x + 1e-15 * I).I * minibatch_x.T * (minibatch_y - p)
            elif self.method == "Gradient":
                accu_gradient = accu_gradient * self.momentum + effective_eta * (
                        -minibatch_x.T * (minibatch_y - p) + self.regularization * self.beta)
                self.beta = self.beta - accu_gradient
            else:
                raise Exception("optimizaation method {:s}".format(self.method))

            loss.append(float(-minibatch_y.T * np.log(p) - (1 - minibatch_y).T * np.log(1 - p)) / len(minibatch_y))

            iterations += 1
            # decay the learning rate
            effective_eta = self.eta / np.power(iterations, self.decay_power)

            if iterations % train_set_size == 0:
                epochs += 1
                mean_loss = np.mean(loss)
                loss = []

                if self.verbose:
                    print("epoch: {:d}. mean loss: {:.6f}. learning rate: {:.8f}".format(epochs, mean_loss, effective_eta))

                if epochs >= self.max_epochs or mean_loss < self.threshold:
                    break

    def predict(self, x):
        x = np.mat(np.c_[[1.0] * x.shape[0], x])
        linear_combination = -np.matmul(x, self.beta)
        linear_combination[linear_combination > 1e2] = 1e2  # prevent overflow
        p = 1.0 / (1.0 + np.power(np.e, linear_combination))

        return p.A1

    def intercept(self):
        return self.beta.A1[0]

    def weights(self):
        return self.beta.A1[1:]

    def batch_loss(self, x, y):
        return

测试该实现在三个数据集,四种优化算法下的表现:

图 4.1 模型在 3 个数据集上的表现

生成上图的测试脚本如下:

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.datasets import make_moons, make_circles, make_classification
from mentat.classification_model import LogisticRegression
from mpl_toolkits.mplot3d import Axes3D

h = .02

names = ["GD No Penalty", "GD L2 Penalty", "Newton", "IRLS"]

classifiers = [
    LogisticRegression(method="Gradient", max_epochs=500, eta=0.01, regularization=0, minibatch_size=100, momentum=0.9, decay_power=0.5, verbose=True),
    LogisticRegression(method="Gradient", max_epochs=500, eta=0.01, regularization=0.01,  minibatch_size=100, momentum=0.9, decay_power=0.5, verbose=True),
    LogisticRegress金融信托配资|配资招商ion(method="Newton", max_epochs=500, eta=0.01,  minibatch_size=100, momentum=0.9, decay_power=0.5, verbose=True),
    LogisticRegression(method="IRLS", max_epochs=500, eta=0.01,  minibatch_size=100, momentum=0.9, decay_power=0.5, verbose=True),
]

X, y = make_classification(n_features=2, n_redundant=0, n_informative=2,
                           random_state=1, n_clusters_per_class=1)
rng = np.random.RandomState(2)
X += 2 * rng.uniform(size=X.shape)
linearly_separable = (X, y)

datasets = [make_moons(noise=0.3, random_state=0), make_circles(noise=0.2, factor=0.5, random_state=1),
            linearly_separable]

figure = plt.figure(figsize=(15, 10))
i = 1

for ds_cnt, ds in enumerate(datasets):

    X, y = ds
    X = StandardScaler().fit_transform(X)
    X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=.4, random_state=42)

    x_min, x_max = X[:, 0].min() - .5, X[:, 0].max() + .5
    y_min, y_max = X[:, 1].min() - .5, X[:, 1].max() + .5
    xx, yy = np.meshgrid(np.arange(x_min, x_max, h),
                         np.arange(y_min, y_max, h))

    cm = plt.cm.RdBu
    cm_bright = ListedColormap(['#FF0000', '#0000FF'])
    ax = plt.subplot(len(datasets), len(classifiers) + 1, i)
    if ds_cnt == 0:
        ax.set_title("Input data")
    ax.scatter(X_train[:, 0], X_train[:, 1], c=y_train, cmap=cm_bright, edgecolors='k')
    ax.scatter(X_test[:, 0], X_test[:, 1], c=y_test, cmap=cm_bright, edgecolors='k', marker="^")
    ax.set_xlim(xx.min(), xx.max())
    ax.set_ylim(yy.min(), yy.max())
    ax.grid(True)
    i += 1

    for name, clf in zip(names, classifiers):
        ax = plt.subplot(len(datasets), len(classifiers) + 1, i,  projection="3d")
        clf.fit(X_train, y_train)
        probability = clf.predict(np.c_[xx.ravel(), yy.ravel()]).reshape(xx.shape)
        ax.plot_surface(xx, yy, probability, rstride=1, cstride=1, alpha=0.6, cmap=cm)
        ax.contourf(xx, yy, probability, zdir='z', offset=-1.0, alpha=0.6, cmap=cm)
        ax.scatter(X_train[:, 0], X_train[:, 1],  -1.0, c=y_train, cmap=cm_bright, edgecolors='k')
        ax.scatter(X_test[:, 0], X_test[:, 1], -1.0, c=y_test, cmap=cm_bright, edgecolors='k', alpha=0.5)
        ax.set_xlim(xx.min(), xx.max())
        ax.set_ylim(yy.min(), yy.max())
        ax.set_zlim(-1.0, 2.0)
        ax.grid(True)
        if ds_cnt == 0:
            ax.set_title(name)
        i += 1

plt.tight_layout()
plt.savefig("lr_contour.png")

做一个广告,作者的基于 numpy 和 pandas 的 python 机器学习库,名叫 Mentat 。欢迎大家交流。

zhangjuefei/mentat​github.com


五、参考书籍

统计学习基础(第2版)(英文) (豆瓣)​book.douban.com最优化导论 (豆瓣)​book.douban.com

关于本文

相关文章


×