逻辑回归算法推理与实现

Overview

逻辑回归通常用于分类算法，例如预测某事是 true 还是 false（二元分类）。例如，对电子邮件进行分类，该算法将使用电子邮件中的单词作为特征，并据此预测电子邮件是否为垃圾邮件。用数学来讲就是指，假设因变量是 Y，而自变量集是 X，那么逻辑回归将预测因变量 \(P(Y=1)\) 作为自变量集 X 的函数。

逻辑回归性能在线性分类中是最好的，其核心为基于样本属于某个类别的概率。这里的概率必须是连续的并且在 (0, 1) 之间（有界）。它依赖于阈值函数来做出称为 Sigmoid 或 Logistic 函数决定的。

学好逻辑回归，需要了解逻辑回归的概念、优势比 (OR) 、Logit 函数、Sigmoid 函数、 Logistic 函数及交叉熵或Log Loss

Prerequisite

odds ratio

explain

odds ratio是预测变量的影响。优势比取决于预测变量是分类变量还是连续变量。

• 连续预测变量：\(OR > 1\) 表示，随着预测变量的增加，事件发生的可能性增加。\(OR < 1\) 表示随着预测变量的增加，事件发生的可能性较小。
• 分类预测变量：事件发生在预测变量的 2 个不同级别的几率；如 A,B，\(OR > 1\) 表示事件在 A 级别的可能性更大。\(OR<1\) 表示事件更低的可能是在A。

例如，假设 X 是受影响的概率，Y 是不受影响的概率，则 \(OR= \frac{X}{Y}\) ，那么 \(OR = \frac{P}{(1-P)}\) ，P是事件的概率。

让概率的范围为 [0,1] ，假设 \(P(success)=0.8\) ，\(Q(failure) = 0.2\) ；\(OR\) 则是 成功概率和失败概率的比值，如：\(O(success)=\frac{P}{Q} = \frac{0.8}{0.2} = 4\) , \(O(failure)=\frac{Q}{P} = \frac{0.2}{0.8} = 0.25\) 。

odds和probability 的区别

• probability 表示在多次实验中，看到改事件的几率，位于 [0,1] 之间

• odds 表示 \(\frac{(事件发生的概率)}{(事件不会发生的概率)}\) 的比率，位于 [0,∞]

例如赛马，一匹马跑 100 场比赛，赢了 80 场，那么获胜的概率是 \(\frac{80}{100} = 0.80 = 80\%\) ，获胜的几率是 \(\frac{80}{20}=4:1\)

总结：probability 和 odds 之间的主要区别：

• “odds”用于描述是否有可能发生事件。相反，probability决定了事件发生的可能性，即事件发生的频率。
• odds以比例表示，probability以百分比形式或小数表示。
• odds通常从 0 ~ ∞ ，其中0定义事件发生的可能性，∞ 表示发生的可能性。相反，probability 介于 0~1之间。因此，probability越接近于0，不发生的可能性就越大，越接近于1，发生的可能性就越高。

通过示例陈述公式

假设一个体校的录取率中，10 个男生中有 7 个被录取，而10 个女生中有3个被录取。找出男生被录取的概率？

那么通过已知条件，设 P 为录取概率，Q则为未被录取的概率，那么

• 男生被录取的概率为：
  • \(P=\frac{7}{10} = 0.7\)
  • \(Q=1-0.7 = 0.3\)
• 女生被录取的概率为：
  • \(P=\frac{3}{10}=0.3\)
  • \(Q=1-0.3=0.7\)
• 录取优势比：
  • \(OR(boy)=\frac{0.7}{0.3}=2.33\)
  • \(OR(Gril) = \frac{0.3}{0.7}=0.42\)

因此，一个男生被录取的几率为 \(OR=\frac{2.33}{0.42}=5.44\)

Logit 函数

logit函数是Odd Ratio 的对数 logarithm , 给出 0~1 范围内的输入，然后将它们转换为整个实数范围内的值。如：假设P，则 \(\frac{P}{(1-P)}\) 为对应的OR；OR 的 logit 的公式为：\(loggit(P) = log(odds) = log(\frac{P}{1-P})\).

以一辆汽车是否出售为例，1为出售，0为不出售，则等式 \(P_i=B_0+B_1 * (Price_i) + \epsilon\)

\(ln(\frac{P}{1-P}) = \beta_0 + \beta_1X_1+\beta_2X_2... + \beta_nX_N\) ,对于简单的逻辑回归，有两个系数：

• \(\beta_0\) 截距 ：X 变量为 0 时的对数 odds ratio
• \(\beta_1\) 斜率：odds ratio随X增加（或减少），1的变化

例如：假设简单逻辑回归模型是 \(Ln(odds) = -5.5 + 1.2*X\) ,那么 \(\beta_0=-5.5\) ，\(\beta_1 = 1.2\) ，意味着，X=0时，\(odds\ ratio = 0\) ，X每增加一个单位 odds ration 增加 1.2（（X 增加2个单位odds ratio增加 2.4....）

求解

通过上面的公式实际上不明白这些具体是什么，就可以通过求P来找到有结果的概率与截距 \(β_0\) 之间的关系，已知 \(n=log_ab\) , $ a^n=b$ ，那么一个简单的逻辑回归公式为 \(log(\frac{P}{1-P}) = \beta_0+\beta1X\) ，对这个公式进行推导：

• \(\frac{P}{1-P} = e^{\beta_0+e^\beta1*X}\)
• \(P = e^{\beta_0+e^\beta1*X} - Pe^{\beta_0+e^\beta1*X}\)
• \(P(1+e^{\beta_0+e^\beta1*X}) = e^{\beta_0+e^\beta1*X}\)
• \(P=\frac{e^{\beta_0+e^\beta1*X}}{1+e^{\beta_0+e^\beta1*X}}\)

当 \(X=0\) ,则 \(\beta_1*X\) 没意义，公式为：\(P = frac{e^{β_0}}{(1+e^{β_0})}\) ，其中e是一个常数，python为 math.e

如果单纯不算概率，只看截距符号，那么满足：

• 如果截距为负号：则产生结果的概率将 < 0.5。
• 如果截距为正号：那么产生结果的概率将 > 0.5。
• 如果截距等于 0：那么得到结果的概率正好是 0.5。

通过例子来说明这点：假设研究为抽烟对心脏健康的影响，下表显示了一个逻辑回归

由表可知，截距为 -1.93，假设smoking系数为0，那么概率带入公式为：\(P=\frac{e^{\beta_0}}{1+e^{\beta_0}} = P=\frac{e^{-1.93}}{1+e^{-1.93}} = 0.126\)(math.e ** -1.93)/(1+math.e ** -1.93)

如果 Smoking是一个连续变量（每年的吸烟量），在这种情况下，Smoking=0 意味着每年使用0公斤烟草的人即不抽烟的人群；那么这个结果就为，不抽烟的人群在未来10年内心脏有问题几率为 0.126。

再如果是吸烟者应该怎么计算，假设，每年吸烟量为3kg，那么公式为：\(P = \frac{e^{β0 + β_1X}}{(1+e^{β0 + β_1X})}\) ，在这里 X=3，那么 \(P=\frac{e^{\beta_1+\beta_2X}}{(1-e^{\beta_1+\beta_2X})} = \frac{e^{-1.93+0.38*3}}{(1-e^{-1.93+0.38*3})} = 0.31\) ；即得出，每年3KG烟草消耗量10年后有心脏问题的概率是 31%

sigmoid

logit 函数的逆函数称Sigmoid 函数，sigmoid方程来源于 logit 为：\(P=\frac{e^{log(odds)}}{(1-e^{log(odds)})} = \frac{1}{e^{-log(odds)+1}} = \frac{1}{1+e^{-z}}\) 。

在python中，np.exp 是求 是求 \(e^{x}\) 的值的函数。正好可以用在sigmod函数中，那么sigmoid可以写为

def sigmoid(z):
    return 1 / (1 + np.exp(-z))

交叉熵或对数损失

交叉熵 Cross-Entropy，通常用于量化两个概率分布之间的差异。用于逻辑回归，公式为：\(H=\sum^{x=n}(P(x) \times log(q(x))\)

Maximum Likelihood Estimation

最大似然估计，Maximum Likelihood Estimation MLE，是概率估算的一种解决方案。MLE在其中寻找一组参数，这些参数将影响数据样本 X 的联合概率的最佳拟合。

首先，定义一个称为 \(\theta\) theta 的参数，该参数定义概率密度函数的选择和该分布的参数。它可能是一个数值向量，其值平滑变化并映射到不同的概率分布及其参数。在最大似然估计中，我们希望在给定特定概率分布及其参数的最大化情况下从联合概率分布中观察数据的概率，形式上表示为：\(P(X|\theta)\) ，在这种情况下，条件概率通常使用分号 ; 而不是竖线 | ，因为 \(\theta\) 不是随机变量，而是未知参数。表达为 \(P(X;\theta)\) ,或 \(P(x_1,x_2,\ ...\ x_n;\theta)\) 。

这样产生的条件概率被称为在给定模型参数 （\(\theta\)）的情况下观察变量 \(X\) 的概率，并使用符号 L 来 表示似然函数。例如：\(L(X;\theta)\)。而最大似然估计的目标是找到使似然函数最大化的一组参数 ( \(\theta\) )，例如产生最大似然值，如：\(max(L(X;\theta))\)

鉴于上述提到的变量 \(X\) 是由n个样本组成，可以将其定义为在给定概率分布参数 \(\theta\) 的情况下，变量 \(X\) 的联合概率,如这里数据样本为 \(x_1,x_2,\ ...\ ,x_n\) 的联合概率，同时表示为 \(L(x_1,x2,\ ...\ ,x_n;\theta)\)

大多数情况下，求解似然方程很复杂。会使用对数似然作为一种解决方案。由于对数函数是单调递增的，因此对数似然和似然中的最优参数是相同的。因此定义条件最大似然估计为：\(log(P(x_i ; h))\)。

用逻辑回归模型替换h，需要假设一个概率分布。在逻辑回归的情况下，假设数据样本为二项式概率分布，其中每个示例都是二项式的一个结果。伯努利分布只有一个参数：成功结果的概率 P，那么为：

• \(P(y=1)=P\)
• \(P(y=0)=1-P\)

那么这个平均值为：\(P(y=1)*1+P(y=0)*0\)，给出P的值公式可以转换为：\(P*1+(1-p)*0\)；这种公式看似没有意义，那么通过一个小例子来了解下

# 二项式似然函数

def likelihood(y, p):
	return p * y + (1 - p) * (1 - y)

# test for y=1
y, p = 1, 0.9
print('y=%.1f, p=%.1f, likelihood: %.3f' % (y, p, likelihood(y, p)))
y, yhat = 1, 0.1
print('y=%.1f, p=%.1f, likelihood: %.3f' % (y, p, likelihood(y, p)))
# test for y=0
y, yhat = 0, 0.1
print('y=%.1f, p=%.1f, likelihood: %.3f' % (y, p, likelihood(y, p)))
y, yhat = 0, 0.9
print('y=%.1f, p=%.1f, likelihood: %.3f' % (y, p, likelihood(y, p)))

# y=1.0, p=0.9, likelihood: 0.900
# y=1.0, p=0.9, likelihood: 0.900
# y=0.0, p=0.9, likelihood: 0.100
# y=0.0, p=0.9, likelihood: 0.100

运行示例会为每个案例打印类别y 和预测概率p，其中每个案例的概率是否接近；这里也可以使用对数更新似然函数，\(log(p) * y + log(1 – p) * (1 – y)\)；最后可以根据数据集中实例求最大似然和最小似然

• \(\sum^{i=1}_n log(p_i) * y_i + log(1 – p_i) * (1 – y_i)\)
• 最小似然使用反转函数，使负对数自然作为最小似然。上面的公式前加 -

对于计算二项式分布的对数似然相当于计算二项式分布[交叉熵，其中P(class)表示第 class 项概率，q() 表示概率分布，\(-(log(q(class0)) \times P(class0) + log(q(class1)) * P(class1))\)

LR算法实例

在研究如何从数据中估计模型的参数之前，我们需要了解逻辑回归准确计算的内容。

模型的线性部分（输入的加权和）计算成功事件的log-odds。

odds ratio：\(\beta_0+\beta_1 \times x_1 + \beta_2 \times x_2\ ...\ \beta_n \times x_n\) 该模型估计了每个级别的输入变量的log-odds。

由上面信息了解到，几率 probability 是输赢的比率 如 1:10 ；probability 可以转换为 odds ratio 即成功概率除以不成功概率：\(or=\frac{P}{1-P}\) ；计算or的对数，被称为log-odds是一种度量单位：\(log(\frac{P}{1-P})\)，而所求的即为 log-odds的逆函数，而在python中 log 函数是对数，求log的逆方法即 exp 返回n的x次方就是log的逆函数。

到这里已经和逻辑回归模型很接近了，对数函数公式可以简化为，\(P=\frac{e^{log(odds)}}{(1-e^{log(odds)})}\) ，以上阐述了如何从log-odds转化为odds，然后在到逻辑回归模型。下面通过Python 中的示例来具体计算 probability 、odds 和 log-odds 之间的转换。假设将成功概率定义为 80% 或 0.8，然后将其转换为odds，然后再次转换为概率。

from math import log
from math import exp

prob = 0.8
print('Probability %.1f' % prob)
# 将 probability 转换为 odds
odds = prob / (1 - prob)
print('Odds %.1f' % odds)
# 将 odds 转换为 log-odds
logodds = log(odds)
print('Log-Odds %.1f' % logodds)
# 转换 log-odds 为  probability
prob = 1 / (1 + exp(-logodds))
print('Probability %.1f' % prob)

# Probability 0.8
# Odds 4.0
# Log-Odds 1.4
# Probability 0.8

通过这个例子，可以看到odds被转换成大约 1.4 的log-odds，然后正确地转回 0.8 的成功概率。

逻辑回归实现

首先将实现分为3个步骤：

• 预测
• 评估系数
• 真实数据集预测

预测

编写一个预测函数，在评估随机梯度下降中的候选系数值时以及在模型最终确定测试数据或新数据进行预测时。

下面是预测predict()函数，它预测给定一组系数的行的输出值。第一个系数是截距，也称为偏差或 b0，它是独立的，不负责输入值。

def predict(row, coefficients):
	p = coefficients[0]
	for i in range(len(row)-1):
		yhat += coefficients[i + 1] * row[i]
	return 1.0 / (1.0 + exp(-p))

准备一些测试数据，Y代表真实的类别

X1					X2						Y
2.7810836 	2.550537003		0
1.465489372	2.362125076		0
3.396561688	4.400293529		0
1.38807019	1.850220317		0
3.06407232	3.005305973		0
7.627531214	2.759262235		1
5.332441248	2.088626775		1
6.922596716	1.77106367		1
8.675418651	-0.242068655	1
7.673756466	3.508563011		1

这里有两个输入值，和三个系数，系数是自定义的固定值，那么预测的公式就为

# 系数为
coef = [-0.406605464, 0.852573316, -1.104746259]
y = 1.0 / (1.0 + e^(-(b0 + b1 * X1 + b2 * X2)))
# 套入公式（sigma）
y = 1.0 / (1.0 + e^(-(-0.406605464 + 0.852573316 * X1 + -1.104746259 * X2)))

完整的代码

# Make a prediction
from math import exp

# Make a prediction with coefficients
def predict(row, coefficients):
	yhat = coefficients[0]
	for i in range(len(row)-1):
		yhat += coefficients[i + 1] * row[i]
	return 1.0 / (1.0 + exp(-yhat))

# test predictions
dataset = [[2.7810836,2.550537003,0],
	[1.465489372,2.362125076,0],
	[3.396561688,4.400293529,0],
	[1.38807019,1.850220317,0],
	[3.06407232,3.005305973,0],
	[7.627531214,2.759262235,1],
	[5.332441248,2.088626775,1],
	[6.922596716,1.77106367,1],
	[8.675418651,-0.242068655,1],
	[7.673756466,3.508563011,1]]
coef = [-0.406605464, 0.852573316, -1.104746259]
for row in dataset:
	yhat = predict(row, coef)
	print("Expected=%.3f, Predicted=%.3f [%d]" % (row[-1], yhat, round(yhat)))

估计系数

这里可以使用我随机梯度下降来估计训练数据的系数值。随机梯度下降需要两个参数：

• 学习率 Learning rate：用于限制每个系数每次更新时的修正量。
• Epochs：更新系数时遍历训练数据的次数。

在每个epoch更新训练数据中每一行的每个系数。系数会根据模型产生的错误进行更新，误差为预期输出与预测值之间的差异。错误会随着epoch增加而减少

将每个都加权，并且这些系数以一致的方式进行更新，用公式可以表示为

b1(t+1) = b1(t) + learning_rate * (y(t) - p(t)) * p(t) * (1 - p(t)) * x1(t)

那么整合一起为

from math import exp

# 预测函数
def predict(row, coefficients):
    p = coefficients[0]
    for i in range(len(row)-1):
        p += coefficients[i + 1] * row[i]
    return 1.0 / (1.0 + exp(-p))

def coefficients_sgd(train, l_rate, n_epoch):
    coef = [0.0 for i in range(len(train[0]))] # 初始一个系数，第一次为都为0
    for epoch in range(n_epoch):
        sum_error = 0
        for row in train:
            p = predict(row, coef)
            # 错误为预期值与实际值直接差异
            error = row[-1] - p
            sum_error += error**2
            # 截距没有输入变量x，这里为row[0]
            coef[0] = coef[0] + l_rate * error * p * (1.0 - p)
            for i in range(len(row)-1):
                # 其他系数更新
                coef[i + 1] = coef[i + 1] + l_rate * error * p * (1.0 - p) * row[i]
        print('>epoch=%d, lrate=%.3f, error=%.3f' % (epoch, l_rate, sum_error))
    return coef

# Calculate coefficients
dataset = [
    [2.7810836,2.550537003,0],
    [1.465489372,2.362125076,0],
    [3.396561688,4.400293529,0],
    [1.38807019,1.850220317,0],
    [3.06407232,3.005305973,0],
    [7.627531214,2.759262235,1],
    [5.332441248,2.088626775,1],
    [6.922596716,1.77106367,1],
    [8.675418651,-0.242068655,1],
    [7.673756466,3.508563011,1]
]
l_rate = 0.3
n_epoch = 100
coef = coefficients_sgd(dataset, l_rate, n_epoch)
print(coef)

# >epoch=92, lrate=0.300, error=0.024
# >epoch=93, lrate=0.300, error=0.024
# >epoch=94, lrate=0.300, error=0.024
# >epoch=95, lrate=0.300, error=0.023
# >epoch=96, lrate=0.300, error=0.023
# >epoch=97, lrate=0.300, error=0.023
# >epoch=98, lrate=0.300, error=0.023
# >epoch=99, lrate=0.300, error=0.022
#[-0.8596443546618897, 1.5223825112460005, -2.218700210565016]

这里跟踪了跟踪每个epoch误差平方的总和，以便我们可以在每个epoch中打印出error，实例中使用 0.3 学习率并训练100 个 epoch，每个epoch会打印出其误差平方，最终会打印总系数集

套用真实数据集

糖尿病数据集 是根据基本的医疗信息，预测印第安人5年内患糖尿病的情况。这是一个二元分类，阴性0与阳性1直接的关系。采用了二项式分布，也可以采用其他分布，如高斯等。

from random import seed
from random import randrange
from csv import reader
from math import exp

# Load a CSV file
def load_csv(filename):
	dataset = list()
	with open(filename, 'r') as file:
		csv_reader = reader(file)
		for row in csv_reader:
			if not row:
				continue
			dataset.append(row)
	return dataset

# Convert string column to float
def str_column_to_float(dataset, column):
	for row in dataset:
		row[column] = float(row[column].strip())

# 找到最小和最大的
def dataset_minmax(dataset):
	minmax = list()
	for i in range(len(dataset[0])):
		col_values = [row[i] for row in dataset]
		value_min = min(col_values)
		value_max = max(col_values)
		minmax.append([value_min, value_max])
	return minmax

# 归一化
def normalize_dataset(dataset, minmax):
	for row in dataset:
		for i in range(len(row)):
			row[i] = (row[i] - minmax[i][0]) / (minmax[i][1] - minmax[i][0])
# k-folds CV实现
def cross_validation_split(dataset, n_folds):
	dataset_split = list()
	dataset_copy = list(dataset)
	fold_size = int(len(dataset) / n_folds)
	for i in range(n_folds):
		fold = list()
		while len(fold) < fold_size:
			index = randrange(len(dataset_copy))
			fold.append(dataset_copy.pop(index))
		dataset_split.append(fold)
	return dataset_split

# 计算准确度百分比
def accuracy_metric(actual, predicted):
	correct = 0
	for i in range(len(actual)):
		if actual[i] == predicted[i]:
			correct += 1
	return correct / float(len(actual)) * 100.0

# 使用CV评估算法
def evaluate_algorithm(dataset, algorithm, n_folds, *args):
	folds = cross_validation_split(dataset, n_folds)
	scores = list()
	for fold in folds:
		train_set = list(folds)
		train_set.remove(fold)
		train_set = sum(train_set, [])
		test_set = list()
		for row in fold:
			row_copy = list(row)
			test_set.append(row_copy)
			row_copy[-1] = None
		predicted = algorithm(train_set, test_set, *args)
		actual = [row[-1] for row in fold]
		accuracy = accuracy_metric(actual, predicted)
		scores.append(accuracy)
	return scores

# 使用系数进行预测
def predict(row, coefficients):
	yhat = coefficients[0]
	for i in range(len(row)-1):
		yhat += coefficients[i + 1] * row[i]
	return 1.0 / (1.0 + exp(-yhat))

# 系数生成
def coefficients_sgd(self, train, l_rate, n_epoch):
    """
    生成系数
    :param train: list, 数据集，可以是训练集
    :param l_rate: float, 学习率
    :param n_epoch:int，epoch，这里代表进行多少次迭代
    :return: None
    """
    coef = [0.0 for i in range(len(train[0]))] # 初始一个系数，第一次为都为0
    for epoch in range(n_epoch):
        sum_error = 0
        for row in train:
            p = self.predict(row, coef)
            # 错误为预期值与实际值直接差异
            error = row[-1] - p
            sum_error += error**2
            # 截距没有输入变量x，这里为row[0]
            coef[0] = coef[0] + l_rate * error * p * (1.0 - p)
            for i in range(len(row)-1):
                # 其他系数更新
                coef[i + 1] = coef[i + 1] + l_rate * error * p * (1.0 - p) * row[i]
                # print('>epoch=%d, lrate=%.3f, error=%.3f' % (epoch, l_rate, sum_error))
            return coef

# 随机梯度下降的逻辑回归算法
def logistic_regression(self, train, test, l_rate, n_epoch):
    predictions = list()
    coef = self.coefficients_sgd(train, l_rate, n_epoch)
    for row in test:
        p = self.predict(row, coef)
        p = round(p)
        predictions.append(p)
    return(predictions)


seed(1)
# 数据预处理
filename = 'pima-indians-diabetes.csv'
dataset = load_csv(filename)
for i in range(len(dataset[0])):
	str_column_to_float(dataset, i)
# 做归一化
minmax = dataset_minmax(dataset)
normalize_dataset(dataset, minmax)
# evaluate algorithm
n_folds = 5
l_rate = 0.1
n_epoch = 100
scores = evaluate_algorithm(dataset, logistic_regression, n_folds, l_rate, n_epoch)
print('Scores: %s' % scores)
print('Mean Accuracy: %.3f%%' % (sum(scores)/float(len(scores))))

# 0.35294117647058826
# Scores: [73.8562091503268, 78.43137254901961, 81.69934640522875, 75.81699346405229, 75.81699346405229]
# Mean Accuracy: 77.124%

上述是对整个数据集的预测百分比，也可以对对应的类的信息进行输出