邏輯迴歸演算法推理與實現

Cylon發表於2022-06-07
邏輯迴歸演算法

Overview

邏輯迴歸通常用於分類演算法，例如預測某事是 true 還是 false（二元分類）。例如，對電子郵件進行分類，該演算法將使用電子郵件中的單詞作為特徵，並據此預測電子郵件是否為垃圾郵件。用數學來講就是指，假設因變數是 Y，而自變數集是 X，那麼邏輯迴歸將預測因變數 \(P(Y=1)\) 作為自變數集 X 的函式。

image

邏輯迴歸效能線上性分類中是最好的，其核心為基於樣本屬於某個類別的概率。這裡的概率必須是連續的並且在 (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，發生的可能性就越高。

Reference

The Difference Between "Probability" and "Odds"

通過示例陳述公式

假設一個體校的錄取率中，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。

通過例子來說明這點：假設研究為抽菸對心臟健康的影響，下表顯示了一個邏輯迴歸

Coefficient Standard Error p-value
Intercept -1.93 0.13 < 0.001
Smoking 0.38 0.17 0.03

由表可知，截距為 -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%

interpret

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}}\) 。

image

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

上述是對整個資料集的預測百分比，也可以對對應的類的資訊進行輸出

Reference

Maximum likelihood estimation

Sigmoid Function

logistic

binary logistic regression

LR implementation

相關文章