机器学习之logistic回归

最新推荐文章于 2022-01-04 15:05:46 发布

董云龙

最新推荐文章于 2022-01-04 15:05:46 发布

阅读量368

点赞数 1

分类专栏：机机器学习文章标签： Logistic回归逻辑回归

本文链接：https://blog.csdn.net/dongyunlon/article/details/79322512

版权

机机器学习专栏收录该内容

4 篇文章 0 订阅

订阅专栏

1. 指数分布族

1.1 定义

指数族分布(exponential family of distributions)亦称指数型分布族，在上世纪30年代中期被提出,在概率学和统计学中,它是一些有着特殊形式的概率分布的集合,是统计中最重要的参数分布族，包含了二项分布、正态分布,指数分布、伯努利分布、泊松分布、gamma分布、beta分布等.指数分布族为很多重要而且常用的概率分布提供了统一的框架,这种一般性有助于表达的方便和从更大的宏观尺度理解这些分布.
指数分布族的概率密度函数PDF或者概率质量函数PMF一般表达式:

P (y; η) = b (y) e x p (η T T (y) - a (η)) (1.1)

$P(y;\eta) = b(y)exp(\eta^TT(y) - a(\eta)) \qquad (1.1)$
分布函数中:
1.

η η $\eta$ 是自然参数(natural parameter), 通常是一个实数,也可以是向量
2.

b(y) b ( y ) $b(y)$ 是底层观测值(underlying measure)
3.

T(y) T ( y ) $T(y)$ 是充分统计量(sufficient statistic),通常为y
4.

a(η) a ( η ) $a(\eta)$ 被称为对数规则化(log normalizer)

当参数a,b,T都固定的时候,就定义了一个以 $\eta$ 为参数的函数族.

$a(\eta)$ 的作用是使PDF的积分等于1,称其为对数规则化,原因如下, 将(1.1)变形:

P (y; η) e x p (a (η)) \cdot P (y; η) = b ( y ) e x p ( η T T ( y ) ) e x p ( a ( η ) ) = b (y) e x p (η T T (y)) (1.2) (1.3)

$\begin{array}{llrr} P( y;\eta ) &= \frac{b(y)exp(\eta^TT(y))}{exp(a(\eta))} &\qquad (1.2) \\ exp(a(\eta))\cdot P(y;\eta) &= b(y)exp(\eta^TT(y)) &\qquad (1.3) \end{array}$
对(1.3)两边同时积分:

e x p (a (η)) \cdot \int P (y; η) = \int b (y) e x p (η T T (y)) (1.4)

$exp(a(\eta))\cdot \int P(y;\eta) = \int b(y)exp(\eta^T T(y)) \qquad (1.4)$
根据定义,概率密度函数PDF的积分为1,所以:

e x p (a (η)) a (η) = \int b (y) e x p (η T T (y)) = ln \int b (y) e x p (η T T (y)) (1.5) (1.6)

$\begin{array}{llrr} exp(a(\eta)) &= \int b(y)exp(\eta^TT(y)) &\qquad (1.5) \\ a(\eta) &= \ln \int b(y)exp(\eta^TT(y)) &\qquad (1.6) \end{array}$
只要

a(η) a ( η ) $a(\eta)$ 满足上述(1.6)中等式,就能保证概率密度函数PDF的积分为1.

1.2 sigmoid函数

假设y满足二项分布,如下:

{P (y = 1) P (y = 0) = θ = 1 - θ (1.7) y \in {0, 1} θ \in [0, 1]

$\left \{ \begin{array}{ll} P(y=1) &= \theta \\ P(y=0) &= 1-\theta \qquad (1.7) \end{array} \right. \\ y \in \{0,1\} \\ \theta \in [0,1]$
将1.7式转化合并为如下形式:

P (y; θ) = θ y \cdot (1 - θ) (1 - y) (1.8)

$P(y;\theta)=\theta^y \cdot (1-\theta)^{(1-y)} \qquad (1.8)$
1.8就是二项分布的概率密度函数,所以,我们可以将其转换为指数分布族的通用框架形式:

P (y; θ) = exp {ln (θ y \cdot (1 - θ) (1 - y)} = exp {y ln θ + (1 - y) ln (1 - θ)} = exp {y ln θ - y ln (1 - θ) + ln (1 - θ)} = exp {y (ln θ - ln (1 - θ)) + ln (1 - θ)} = exp {y (ln θ ( 1 - θ ) ) + ln (1 - θ)}

$\begin{array}{ll} P(y;\theta) &= \exp \{\ln(\theta^y \cdot (1-\theta)^{(1-y)}\} \\ &= \exp \{ y\ln \theta + (1-y)\ln (1-\theta) \} \\ &= \exp\{ y\ln \theta -y\ln (1-\theta) + \ln(1-\theta)\} \\ &= \exp\{ y(\ln \theta -\ln (1-\theta)) + \ln (1-\theta)\} \\ &= \exp\{ y(\ln \frac{\theta} {(1-\theta))} + \ln(1-\theta)\} \\ \end{array}$
和1.1公式比较,可知:

b (y) η T (y) a (η) = 1 = ln θ 1 - θ (1.9) = y = ln (1 - θ) (1.10)

$\begin{array}{ll} b(y) &= 1 \\ \eta &= \ln \frac{\theta}{1-\theta} \qquad (1.9)\\ T(y) &= y \\ a(\eta) &= \ln(1-\theta) \qquad (1.10) \end{array}$
由公式1.9可知:

θ = 1 1 + e - η (1.11)

$\theta = \frac{1}{1+e^{-\eta}} \qquad (1.11)$
由公式1.11可知,二项分布的概率参数

θ θ $\theta$ `和指数分布族的自然参数`

η η $\eta$ 之间的关系.
公式1.11就是sigmoid函数. 将式1.11带入式1.10,可知:

a (η) = ln (1 + e η)

$a(\eta) = \ln (1+ e^\eta)$

2. Logstic回归

2.1 推导模型

假设:
1. 样本集和样本:

X x (i) = [x (1), x (2), x (3), \dots, x (m)] = [x (i) 1, x (i) 2, x (i) 3, \dots, x (i) n] T

$\begin{array}{ll} X &= [x^{(1)}, x^{(2)}, x^{(3)}, \cdots, x^{(m)}] \\ x^{(i)} &= [x^{(i)}_1, x^{(i)}_2, x^{(i)}_3, \cdots, x^{(i)}_n]^T \end{array}$
2. 假设分类y在x给定的情况下,满足二项分布,并且其为1的概率是样本特征x的函数,参数

θ θ $\theta$ :

P (y = 1 | x; θ) P (y = 0 | x; θ) P (y | x; θ) y h θ (x) = h θ (x) = 1 - h θ (x) = h θ (x) y \cdot (1 - h θ (x)) 1 - y (2.1) \in {0, 1} \in {0, 1}

$\begin{array}{rl} P(y=1|x;\theta) &= h_\theta(x) \\ P(y=0|x;\theta) &= 1 - h_\theta(x) \\ P(y|x;\theta) &= h_\theta(x)^y \cdot (1 - h_\theta(x))^{1-y} \qquad (2.1)\\ y &\in \{0, 1\} \\ h_\theta(x) &\in \{0, 1\} \end{array}$
已知二项分布满足指数分布族通用框架,在其他参数固定的情况下,不同的自然参数

η η $\eta$ 代表着不同概率的二项分布,所以,我们假设自然参数和样本特征之间存在线性的关系:

η η = θ 1 x 1 + θ 2 x 2 + \dots + θ n x n = θ T x (2.2)

$\begin{array}{ll} \eta &= \theta_1x_1 + \theta_2x_2 + \cdots +\theta_nx_n \\ \eta &= \theta^Tx \qquad (2.2) \end{array}$
结合公式2.5,可知:

h θ (x) = 1 1 + e - θ T x (2.3)

$h_\theta(x) = \frac{1}{1+e^{-\theta^Tx}} \qquad (2.3)$
由于公式2.2的假设, 所以,Logstic回归也属于广义线性模型, 这里并不直接假设y和x之间的关系,否则就是线性模型,而是假设在已知x的情况下,y的取值存在一个概率,x和这个概率之间存在关系,而

η η $\eta$ `又代表了不同概率的同一分布,所以很自然假设x和`

η η $\eta$ 之间存在线性关系,并由此来训练模型.

2.2 构建模型

已知y的概率分布,将其应用到所有样本,求取参数向量 $\theta^T$ 的值,使其应当满足,所有样本的概率之积最大,即假设,已经收集到的样本是最有可能出现的, 未收集到的样本出现的概率相对较低.

L (θ) = P (y (1) | x (1); θ) \cdot P (y (2) | x (2); θ) \dots P (y (m) | x (m); θ) = \prod m i = 0 h θ (x (i)) y (i) \cdot (1 - h θ (x (i))) 1 - y (i) (2.4)

$\begin{array}{ll} L(\theta) &= P(y^{(1)}|x^{(1)};\theta) \cdot P(y^{(2)}|x^{(2)};\theta) \cdots P(y^{(m)}|x^{(m)};\theta) \\ &= \prod_{i=0}^m h_\theta(x^{(i)})^{y^{(i)}} \cdot (1-h_\theta(x^{(i)}))^{1-y^{(i)}} \qquad (2.4) \end{array}$
如上,根据需求,应该求解2.4关于

θ θ $\theta$ 的函数,使其结果最大.由于连乘积较难计算最大值,所以,左右两边同时求对数,因为对数函数递增,不影响原函数的升降.

l (θ) = ln \prod m i = 0 h θ (x (i)) y (i) \cdot (1 - h θ (x (i))) 1 - y (i) = \sum m i = 0 y (i) ln h θ (x (i)) + (1 - y (i)) ln (1 - h θ (x (i))) (2.5)

$\begin{array}{ll} l(\theta) &= \ln \prod_{i=0}^m h_\theta(x^{(i)})^{y^{(i)}} \cdot (1-h_\theta(x^{(i)}))^{1-y^{(i)}} \\ &= \sum_{i=0}^m y^{(i)}\ln h_\theta(x^{(i)}) +(1-y^{(i)})\ln (1-h_\theta(x^{(i)})) \qquad (2.5) \end{array}$
为了求解式2.5的最大值,这里选择使用 梯度上升法, 顾名思义,和梯度下降法相反,不断地叠加梯度,进而不断靠近极值点.第一步需要求

θj θ j $\theta_j$ 的偏导数:

\partial l ( θ j ) \partial θ j = \sum m i = 0 y (i) 1 h θ ( x ( i ) ) \partial h θ ( x ( i ) ) \partial θ j - (1 - y (i)) 1 1 - h θ ( x ( i ) ) \partial h θ ( x ( i ) ) \partial θ j = \sum m i = 0 \partial h θ ( x ( i ) ) \partial θ j (( 1 - h θ ( x ( i ) ) y ( i ) - h θ ( x ( i ) ) ( 1 - y ( i ) ) h θ ( x ( i ) ) ( 1 - h θ ( x ( i ) ) )) = \sum m i = 0 \partial h θ ( x ( i ) ) \partial θ j (y ( i ) - h θ ( x ( i ) ) h θ ( x ( i ) ) ( 1 - h θ ( x ( i ) ) )) (2.6)

$\begin{array}{ll} \frac{\partial l(\theta_j)}{\partial \theta_j} &= \sum_{i=0}^m y^{(i)} \frac{1}{h_\theta(x^{(i)})} \frac{\partial h_\theta(x^{(i)})}{\partial \theta_j} - (1-y^{(i)}) \frac{1}{1-h_\theta(x^{(i)})} \frac{\partial h_\theta(x^{(i)})}{\partial \theta_j} \\ &= \sum_{i=0}^m \frac{\partial h_\theta(x^{(i)})}{\partial \theta_j} (\frac{(1-h_\theta(x^{(i)})y^{(i)} - h_\theta(x^{(i)})(1-y^{(i)})}{h_\theta(x^{(i)}) (1-h_\theta(x^{(i)}))}) \\ &= \sum_{i=0}^m \frac{\partial h_\theta(x^{(i)})}{\partial \theta_j}(\frac{y^{(i)} - h_\theta(x^{(i)})} {h_\theta(x^{(i)}) (1-h_\theta(x^{(i)}))}) \qquad (2.6) \end{array}$
到这里,就需要计算

hθ(x) h θ ( x ) $h_\theta(x)$ `关于`

θj θ j $\theta_j$ 的偏导数了:

\partial h θ ( x ) \partial θ j h θ (x (i)) (1 - h θ (x (i))) = (1 1 + e - θ T x)' θ j = - (1 + e - θ T x) - 2 e - θ T x \cdot (- θ T x)' θ j = (1 + e - θ T x) - 2 e - θ T x \cdot x j (2.7) = (1 1 + e - θ T x) \cdot (1 - 1 1 + e - θ T x) = (1 1 + e - θ T x) \cdot (e - θ T x 1 + e - θ T x) = (1 + e - θ T x) - 2 e - θ T x (2.8)

$\begin{array}{ll} \frac{\partial h_\theta(x)}{\partial \theta_j} &= (\frac{1}{1+ e^ {-\theta^Tx}})'_{\theta_j} \\ &= -(1+e^{-\theta^Tx})^{-2}e^{-\theta^Tx} \cdot(-\theta^Tx)'_{\theta_j} \\ &= (1+e^{-\theta^Tx})^{-2}e^{-\theta^Tx} \cdot x_j \qquad (2.7) \\ h_\theta(x^{(i)}) (1-h_\theta(x^{(i)})) &= (\frac{1}{1+ e^ {-\theta^Tx}}) \cdot (1-\frac{1}{1+ e^ {-\theta^Tx}}) \\ &= (\frac{1}{1+ e^ {-\theta^Tx}}) \cdot (\frac{e^ {-\theta^Tx}}{1+ e^ {-\theta^Tx}}) \\ &=(1+e^{-\theta^Tx})^{-2}e^{-\theta^Tx} \qquad (2.8) \end{array}$
将式2.7和式2.8带入公式2.6中,可得:

\partial l ( θ j ) \partial θ j = \sum i = 0 m (y (i) - h θ (x (i))) \cdot x (i) (2.9)

$\frac{\partial l(\theta_j)}{\partial \theta_j} = \sum_{i=0}^m (y^{(i)} - h_\theta(x^{(i)}))\cdot x^{(i)} \qquad (2.9)$
结果很优美,就是分类y减去预测分类

hθ(x) h θ ( x ) $h_\theta(x)$ 的差值,乘以x.
至此,我们能够知道参数

θj θ j $\theta_j$ 的更新规则:

θ j = θ j + α \sum i = 0 m (y (i) - h θ (x (i))) \cdot x (i); j \in {1, 2, \dots, n}, α > 0 (2.10)

$\theta_j = \theta_j + \alpha \sum_{i=0}^m(y^{(i)} - h_\theta(x^{(i)}))\cdot x^{(i)} ;\qquad j \in \{1, 2, \cdots, n\}, \quad \alpha >0 \qquad (2.10)$
其中,

α α $\alpha$ `是步长参数,用以控制每次更新`

θ θ $\theta$ 的增量大小,也叫学习速率,这个需要在实践中进行调整,其值过小会导致多次迭代才能收缩,过大会导致越过最优点发生震荡现象.
梯度下降会导致局部极值点的产生,解决方法是随机初始化,寻找多个最优点结果,在这些最优解中寻找最终结果.但是,对于本例,再次对公式2.9求导可知,二次导数小于0,属于凸函数(看定义,有些教材是二次导数大于0,称为凸函数,凸凹函数一样,只有一个最优解),最优解只有一个.
当我们使用Logistic回归时,只需要使用公式2.3和公式2.10即可.

3. 代码实现

3.1 伪代码

梯度下降法(梯度上升法)有两种实现方式:
1. 批梯度下降法 (batch gradient descent)

更新每一参数时,都使用所有的样本进行更新.

Repeat Until Converge{

θ j = θ j + α \sum i = 0 m (y (i) - h θ (x (i))) \cdot x (i) (f o r e v e r y j \in {1, 2, \dots, n})

$\theta_j = \theta_j + \alpha \sum_{i=0}^m(y^{(i)} - h_\theta(x^{(i)}))\cdot x^{(i)} \\ (for \quad every \quad j \in \{1, 2, \cdots, n\})$
}

优点:全局最优解;迭代次数相对较少;易于并行实现;
缺点:当样本数目很多时,训练过程会很慢
2. 随机梯度下降法(stochastic gradient descent) 或者增量梯度下降法(increment gradient descent)

为了解决样本数据量过大的问题,提出了随机梯度下降法,通过每个样本来迭代更新一次,不需要每次都计算全部样本:

Repeat Until Converge{
for i = 1,2,…m{

θ j = θ j + α (y (i) - h θ (x (i))) \cdot x (i) (f o r e v e r y j \in {1, 2, \dots, n})

$\theta_j = \theta_j + \alpha (y^{(i)} - h_\theta(x^{(i)}))\cdot x^{(i)} \\ (for \quad every \quad j \in \{1, 2, \cdots, n\})$
}
}

如果样本量很大的情况下,例如十几万,可能只需要其中的几万条或者几千条,就已经迭代到最优解,对比上面的批量梯度下降，迭代一次需要用到十几万训练样本，一次迭代不可能最优，如果迭代10次的话就需要遍历训练样本10次。但是，SGD伴随的一个问题是噪音较BGD要多，使得SGD并不是每次迭代都向着整体最优化方向, 如果连续出现相似的样本,导致迭代结果变化不大,迭代可能就此终止在该局部最优解处.
优点:训练速度快
缺点:准确度下降，并不是全局最优；不易于并行实现.
既然如此,实际上还有另外一种方案,结合1,2方案,每次迭代选取10个样本,一直迭代至收缩,这种方法称为小批量梯度下降法(MBGD),即保证训练速度比较快,又保证最终参数的准确率,理解1,2方案,该方法实现就很简单了,不详述.

3.2 python实现

# -- coding: utf-8 --
import numpy as np
import matplotlib.pyplot as plt


def func_sigmoid(x):
    sig = 1.0 / (1 + np.exp(-x))
    return sig


def draw_split_line(theta):
    """
    hθ(z) = logistic(z)
    z = θT * x  z:特征空间  x:输入空间
    θT * X =0 唯一确定一个超平面L, 使数据线性可分, 此时 L平面上方数据对应 z > 0, L平面下方数据对应z < 0
    :param theta:
    :return:
    """
    x = np.linspace(-10, 10, 50)
    vector_x = np.row_stack((np.ones((1, 50)), np.mat(x)))
    y = -((np.dot(theta[0, 0:2], vector_x)) / theta[0, 2])
    plt.plot(x, y.T)


def draw_train_data(train_x, train_y):
    num_futures, num_samples = np.shape(train_x)
    for i in range(0, num_samples):
        color = 'or' if train_y[i, 0] == 1 else 'ob'
        plt.plot(train_x[1, i], train_x[2, i], color, markersize=5)


def get_train_data():
    train_x = []
    train_y = []
    data_file = open("../doc/logistic_train_data.txt")
    for line in data_file.readlines():
        line_array = line.strip().split()
        train_x.append([1.0, float(line_array[0]), float(line_array[1])])
        train_y.append(int(line_array[2]))
    return np.mat(train_x).transpose(), np.mat(train_y).transpose()


def train(theta, train_x, train_y, times, alpha, precision, train_type):
    """
    Parameters
    ----------
    theta: 目标参数
    train_x:训练集 x维度
    train_y:训练集 y维度
    times:最大循环次数
    alpha:步长
    precision:精度,满足精度时停止循环
    train_type: 训练模型, base: 全量梯度下降, incremental:随机梯度/增量梯度下降法
    ----------
    """
    for c in range(0, times):
        theta_bak = theta.copy()
        if train_type == 'base':
            theta = train_mode_base(theta, train_x, train_y, alpha)
        elif train_type == 'incremental':
            theta = train_mode_incremental(theta, train_x, train_y, alpha)
        error = theta_bak - theta
        # 判断精度是否满足要求
        if (np.fabs(error) < precision).any():
            print "precision:", precision, "count:", c, "error:", error
            break
    return theta


# 全量梯度下降法
def train_mode_base(theta, train_x, train_y, alpha):
    num_futures, num_samples = np.shape(train_x)
    for j in range(0, num_futures):
        step_sum = 0
        for i in range(0, num_samples):
            train_xi = train_x[:, i]
            train_xji = train_x[j, i]
            train_yi = train_y[i, 0]
            theta_x = np.dot(theta, train_xi)[0, 0]
            sig = func_sigmoid(theta_x)
            error = train_yi - sig
            step = error * train_xji
            step_sum += step
        theta[0, j] += alpha * step_sum
    return theta


# 增量/随机梯度下降法
def train_mode_incremental(theta, train_x, train_y, alpha):
    num_futures, num_samples = np.shape(train_x)
    for i in range(0, num_samples):
        for j in range(0, num_futures):
            train_xi = train_x[:, i]
            train_xji = train_x[j, i]
            train_yi = train_y[i, 0]
            theta_x = np.dot(theta, train_xi)[0, 0]
            sig = func_sigmoid(theta_x)
            error = train_yi - sig
            step = error * train_xji
            theta[0, j] += alpha * step
    return theta


def main():
    plt.figure("mode")
    train_x, train_y = get_train_data()
    # 绘制样本数据
    draw_train_data(train_x, train_y)
    # 训练
    theta = np.mat([1.0, 1.0, 1.0])
    theta = train(theta, train_x, train_y, 1000, 0.1, 0.001, 'base')
    print theta
    draw_split_line(theta)
    plt.show()

if __name__ == '__main__':
    main()

训练数据下载地址: logistic_train_data.txt

参考:
1. https://www.cnblogs.com/maybe2030/p/5089753.html
2. http://blog.csdn.net/saltriver/article/details/55105285
3. http://blog.csdn.net/stdcoutzyx/article/details/9113681
4. http://blog.csdn.net/stdcoutzyx/article/details/9207047

董云龙

关注

1
点赞
踩
1

收藏

觉得还不错? 一键收藏
0
评论
机器学习之logistic回归

1. 指数分布族1.1 定义指数族分布(exponential family of distributions)亦称指数型分布族，在上世纪30年代中期被提出,在概率学和统计学中,它是一些有着特殊形式的概率分布的集合,是统计中最重要的参数分布族，包含了二项分布、正态分布,指数分布、伯努利分布、泊松分布、gamma分布、beta分布等.指数分布族为很多重要而且常用的概率分布提供了统一的框架...
复制链接

扫一扫