邏輯斯蒂迴歸與最大熵模型初探

1. 演算法概述

0x1：邏輯斯蒂迴歸

邏輯斯蒂迴歸（logistic regression）是統計學習中的經典分類方法，它屬於一種對數線性模型（轉化為對數形式後可轉化為線性模型）

從概率角度看，邏輯斯蒂迴歸本質上是給定特徵條件下的類別條件概率

0x2：最大熵準則

最大熵模型的原則：

承認已知事物（知識）；
對未知事物不做任何假設，沒有任何偏見。

對一個隨機事件的概率分佈進行預測時，我們的預測應當滿足全部已知條件，而對未知的情況不要做任何主觀假設。在這種情況下，概率分佈最均勻，預測的風險最小。

因為這時概率分佈的資訊熵最大，所以人們把這種模型叫做“最大熵模型”（Maximum Entropy）

0x3：邏輯斯蒂迴歸和深度神經網路DNN的關係

邏輯斯蒂迴歸模型中的sigmoid函式輸出是一個概率置信度，我們把單層邏輯斯蒂迴歸的條件概率輸出直接進行argmax判斷就是傳統邏輯斯蒂迴歸模型，如果stacking層層堆疊起來，就成為了一個人工神經網路

Relevant Link:

https://www.cnblogs.com/little-YTMM/p/5582271.html

 

1. 邏輯斯蒂迴歸模型

0x1：邏輯斯蒂分佈（logistic distribution）

設 X 是隨機變數，X 服從邏輯斯蒂分佈是指 X 具有下列分佈函式和密度函式

式中，為位置引數，為形狀引數

邏輯斯蒂迴歸密度函式；和分佈函式圖形如下圖所示

分佈函式屬於邏輯斯蒂函式，其圖形是一條S形曲線（sigmoiid curve），該曲線以點為中心對稱

曲線在中心附近增長變化較快，在兩端增長變化較慢。形狀引數越小，sigmoid曲線在中心附近增長得越快

0x2：二項邏輯斯蒂迴歸模型 - 二分類

二項邏輯斯蒂迴歸模型（binominal logistic regression model）是一種分類模型，由條件概率分佈表示，形式為引數化（概率離散化）的邏輯斯蒂分佈

對於給定的輸入例項 x，按照上式可以得到和。二項邏輯斯蒂迴歸比較兩個條件概率值的大小（類似樸素貝葉斯的最大後驗概率比較），將例項 x 分到概率較大的那一類

1. 藉助二項邏輯斯蒂迴歸來驗證一個論斷：邏輯斯蒂迴歸是對數線性模型

我們來考察一下邏輯斯蒂迴歸模型的特點。一個事件的機率（odds）是指該事件發生的概率與該事件不發生的概率的比值，對二項邏輯斯蒂迴歸而言，機率公式為：

這就是說，在邏輯斯蒂迴歸模型中，輸出 Y = 1的對數機率是輸入 x 的線性函式，或者說，輸出 Y = 1的對數機率是由輸入 x 的線性函式表示的模型，即邏輯斯蒂迴歸模型

值得注意的是，這個論斷對多項情況下也是成立的，只是證明更復雜一些

0x3：多項式邏輯斯蒂迴歸模型 - 多分類

多項式邏輯斯蒂迴歸模型（multi-norminal logistic regression model）用於多分類。設離散型隨機變數 Y 的取值集合是｛1，2，....，K｝，那麼多項邏輯斯蒂迴歸模型是：

 

2. 最大熵模型

0x1： 最大熵原理

最大熵原理是概率模型中的一個準則。最大熵原理認為，學習概率模型時，在所有可能的概率模型（分佈）中，熵最大的模型是最好的模型。通常用約束條件來確定概率模型的集合，所以最大熵原理也可以表述為在滿足約束條件的模型集合中選取熵最大的模型

假設離散隨機變數 X 的概率分佈是 P（X），則其熵為：

熵滿足下列不等式：，式中是 X 的取值個數，當且僅當 X 的分佈是均勻分佈時右邊的等號成立，也即當 X 服從均勻分佈時，熵最大

直觀上說：

最大熵原理認為要選擇的概率模型首先必須滿足已有的事實，即約束條件，在沒有更多資訊的情況下，那些不確定的部分都是"等可能的"。最大熵原理通過熵的最大化來表示等可能性，因為"等可能"不容易工程化，而熵則是一個可優化的數值指標

1. 用一個例子來說明最大熵原理

假設隨機變數 X 有5個取值｛A，B，C，D，E｝，要估計各個取值的概率

解：這些概率滿足以下約束條件：

但問題是，滿足這個約束條件的概率分佈有無窮多個，如果沒有任何其他資訊，要對概率分佈進行估計，一個辦法就是認為這個分佈中的各個取值是等概率的：

等概率表示了對事實的無知，因為沒有更多的資訊，這種判斷是合理的。

但是有時，能從一些先驗知識中得到一些對概率值的約束條件，例如：

滿足這兩個約束條件的概率分佈仍然有無窮多個，在缺失其他資訊的情況下，可以繼續使用最大熵原則，即：

我們在樸素貝葉斯法的最大似然估計中，在沒有明確先驗知識的情況下，對先驗概率設定等概率分佈（即熵值最大）的做法就是符合"最大熵原理"的一種做法。但是在引數估計前如果對先驗分佈有一定的知識，就可以變成MAP估計

0x2：最大熵模型的定義

最大熵原理是統計學習的一般原理，將它應用到分類得到最大熵模型。假設分類模型是一個條件概率分佈，表示給定的輸入Ｘ，以條件概率輸出Ｙ

給定一個訓練資料集：

學習的目標是用最大熵原理選擇最好的分類模型。

我們先來定義約束條件，即：

特徵函式關於經驗分佈的期望值（從樣本中估計），與，特徵函式關於模型與經驗分佈的期望值。如果模型能夠獲取訓練資料中的資訊，那麼就可以假設這兩個期望相等：

我們將上式作為模型學習的約束條件，加入有 n 個特徵函式

，那麼就有 n 個約束條件

1. 最大熵模型數學定義

 假設滿足左右約束條件的模型集合為：

定義在條件概率分佈上的條件熵為：

則模型集合C中條件熵最大的模型成為最大熵模型

 

3. 邏輯斯蒂迴歸模型及最大熵模型策略

0x1：約束最優化問題 - 解決最大熵模型空間的搜尋問題

最大熵模型的學習過程就是求解最大熵模型的過程，最大熵模型的學習可以形式化為約束最優化問題

對於給定的訓練資料集以及特徵函式，最大熵模型的學習等價於約束最優化問題：

按照對偶等價的思想，將求最大值問題改寫為等價的求最小值問題：

求解上式約束最優化問題所得出的解，就是最大熵模型學習的解。下面進行一個推導

將約束最優化問題的原始問題轉換為無約束最優化的對偶問題，通過解對偶問題求解原始問題。首先，引進拉格朗日乘子w0，w1，...，wn，定義拉格朗日函式：L（P，w）

最優化的原始問題是：

對偶問題是：

首先，求解對偶問題內部的極小化問題：，它是 w 的函式，將其記作：，稱為對偶函式。同時，將其解記作：

具體地，求對的偏導數：

令偏導數等於0，在的情況下，解得：

由於，得

，其中：

稱為規範化因子；是特徵函式；是特徵的權值。上二式表示的就是最大熵模型，這裡 w 是最大熵模型中的引數向量

之後，求解對偶問題外部的極大化問題：，將其解記為，即

1. 用一個例子來說明最大熵原理

文章上面我們通過直接應用最大熵原理和直觀上的直覺解決了一個在有限約束下的模型引數估計問題，這裡我們繼續運用拉格朗日乘子法求解約束最優化的思想再次求解該問題

分別以y1，y2，y3，y4，y5表示A，B，C，D，E。於是最大熵模型學習的最優化問題是：

引進拉格朗日乘子w0，w1，定義拉格朗日函式：

根據拉格朗日對偶性，可以通過求解對偶最優化問題得到原始最優化問題的解：

首先求解關於P的極小化問題，為此，固定w1，w1，求偏導數

令各偏導數等於0，解得：

於是：

再求解關於 w 的極大化問題：

分別對w0和w1求偏導數並令其為0，得到：

可以看到，結果和我們直觀上理解的是一致的

0x2：最大似然估計 - 經驗風險最小化擬合訓練樣本

上一小節我們看到，最大熵模型是由

和 表示的條件概率分佈。

下面我們來證明其對偶函式的極大化等價於最大熵模型的極大似然估計（即最大熵原理和極大似然原理的核心思想是互通的）

已知訓練資料的經驗概率分佈，條件概率分佈的對數似然函式表示為：

當條件概率分佈是最大熵模型時，對數似然函式為：

再看條件概率最大熵模型的對偶函式

比較上面兩式可得：

即對偶函式等價於對數似然函式，於是證明了最大熵模型學習中的對偶函式極大化等價於最大熵模型的極大似然估計這個論斷

這樣，最大熵模型的學習問題就轉換為具體求解對數似然函式極大化或對偶函式極大化的問題

可以將最大熵模型寫成更一般的形式：，這裡，為任意實值特徵函式

0x3：邏輯斯蒂迴歸與最大熵模型的等價性

最大熵模型與邏輯斯蒂迴歸模型有類似的形式，它們又稱為對數線性模型（log linear model）。

它們在演算法策略上也存在共通的形式和思想，模型學習都是在給定的訓練資料條件下對模型進行極大似然估計（遵循最大熵原理），或正則化的極大似然估計

在最大熵模型中, 令:

即可得到邏輯斯蒂迴歸模型.

Relevant Link:

http://blog.csdn.net/foryoundsc/article/details/71374893

 

4. 邏輯斯蒂迴歸、最大熵演算法 - 引數估計

我們既然已經有了最大熵原理和最大似然估計了，不就可以直接得到目標函式（或條件概率）的全域性最優解了嗎？但是就和決策樹學習中一樣，最大熵增益原理只是指導思想，具體到工程化實現時還需要有一個高效快捷的方法來求得"最佳引數"

邏輯斯蒂迴歸模型、最大熵模型學習歸結為以似然函式為目標函式的最優化問題，因為目標函式中包含有大量的乘法項直接求導求極值十分困難，因此通過通過迭代演算法求解（這點在深度神經網路中也是一樣）。從最優化的觀點看，這時的目標函式具有很好的性質，它是光滑的凸函式，因此多種最優化的方法都適用，保證能找到全域性最優解，常用的方法有改進的迭代尺度法、梯度下降法、牛頓法、擬牛頓法

0x1：改進的迭代尺度法

改進的迭代尺度法（improved interative scaling IIS）是一種最大熵模型學習的最優化演算法。

已知最大熵模型為：

對數似然函式為：

目標是通過極大似然估計學習模型引數，即求對數似然函式的極大值。因為似然函式 = 概率密度（而這裡的概率密度函式就是最大熵模型），所以求對數似然函式的極大值就是在求最大熵模型的最優解

IIS的想法是：

假設最大熵模型當前的引數向量是，我們希望找到一個新的引數向量，使得模型的對數似然函式值增大。如果能有這樣一種引數向量的更新方法，那麼就可以重複使用這一方法，（逐步迭代），直至找到對數似然函式的最大值

對於給定的（本輪迭代）經驗分佈，模型引數從到，對數似然函式的改變數是：

利用不等式：

建立對數似然函式改變數的下界：

將右端記為：

於是有：

即是對數似然函式本輪迭代改變數的一個下界。如果能找到適當的使下界提高，那麼對數似然函式也會提高。

然而，函式中的是一個向量，含有多個變數，不易同時優化。IIS試圖一次只優化其中一個變數，而固定其他變數

為達到這一目的，IIS進一步降低下界。具體地，IIS引進一個量

因為是二值函式，故表示所有特徵在（x，y）出現的次數，這樣，可以改寫為：

根據Jensen不等式上式可改寫為：

記不等式右端為：

於是得到：

這裡，是對數似然改變數的一個新的（相對不緊）下界。於是求對的偏導數：

上式中，除了外不包含任何其他變數，令偏導數為0得到：

依次對求解方程可以求出

演算法的過程可以參閱下面的python程式碼實現

 

5. 邏輯斯蒂迴歸模型和感知機模型的聯絡與區別討論

我們知道，感知機模型對輸入 x 進行分類的線性函式，其值域（輸入）為整個實屬域，以一元（1維特徵變數）感知機為例，感知機的判別函式實際上是一個階躍函式

邏輯斯蒂迴歸模型相當於改進版本的感知機，它在感知機的輸出後面加上了一個sigmoid函式，使其概率分佈函式增加了一些新的特性

1. f(x;β)∈[0,1].
2. f應該至少是個連續函式. 這是因為我們希望模型f的輸出能夠隨 x平滑地變化.
3. f應該儘可能簡單.

sigmoid不僅完成了值域概率化壓縮，同時還使概率分佈函式連續可導，通過邏輯斯蒂迴歸模型可以將線性函式轉換為概率：

這時，線性函式的值（輸入）越接近正無窮，概率值就越接近於1；線性函式的值越接近於負無窮，概率值就越接近於0，即實現了判別結果概率化壓縮，這樣的模型就是邏輯斯蒂迴歸模型

 

6. 用Python實現最大熵模型（MNIST資料集） 

0x1：資料地址

https://github.com/LittleHann/lihang_book_algorithm/blob/master/data/train_binary.csv

0x2：code

# -*- coding: utf-8 -*-

import pandas as pd
import numpy as np
import os

import time
import math
import random

from collections import defaultdict

from sklearn.cross_validation import train_test_split
from sklearn.metrics import accuracy_score


class MaxEnt(object):

    def init_params(self, X, Y):
        self.X_ = X
        self.Y_ = set()

        # 計算(x, y)和(x)的總數
        self.cal_Pxy_Px(X, Y)

        self.N = len(X)                 # 訓練集大小
        self.n = len(self.Pxy)          # 書中(x,y)總數n
        self.M = 10000.0                # 書91頁那個M，但實際操作中並沒有用那個值
        # 可認為是學習速率

        self.build_dict()

        # 計算特徵函式f(x,y)關於經驗分佈P~(X,Y)的期望值
        self.cal_EPxy()

    def build_dict(self):
        self.id2xy = {}
        self.xy2id = {}

        for i, (x, y) in enumerate(self.Pxy):
            self.id2xy[i] = (x, y)
            self.xy2id[(x, y)] = i

    def cal_Pxy_Px(self, X, Y):
        self.Pxy = defaultdict(int)
        self.Px = defaultdict(int)

        for i in xrange(len(X)):
            x_, y = X[i], Y[i]
            self.Y_.add(y)

            for x in x_:
                self.Pxy[(x, y)] += 1
                self.Px[x] += 1

    def cal_EPxy(self):
        '''
        計算書中82頁最下面那個期望
        '''
        self.EPxy = defaultdict(float)
        for id in xrange(self.n):
            (x, y) = self.id2xy[id]
            # 特徵函式f(x,y)是0/1指示函式，所以頻數/總數 = 概率P(x,y)
            self.EPxy[id] = float(self.Pxy[(x, y)]) / float(self.N)

    def cal_pyx(self, X, y):
        result = 0.0
        for x in X:
            if self.fxy(x, y):
                id = self.xy2id[(x, y)]
                result += self.w[id]
        return (math.exp(result), y)

    def cal_probality(self, X):
        '''
        計算書85頁公式6.22
        '''
        # 計算P(Y|X)條件概率
        Pyxs = [(self.cal_pyx(X, y)) for y in self.Y_]
        Z = sum([prob for prob, y in Pyxs])
        return [(prob / Z, y) for prob, y in Pyxs]

    def cal_EPx(self):
        '''
        計算書83頁最上面那個期望
        '''
        self.EPx = [0.0 for i in xrange(self.n)]

        for i, X in enumerate(self.X_):
            # 得到P(y|x)條件概率
            Pyxs = self.cal_probality(X)

            for x in X:
                for Pyx, y in Pyxs:
                    if self.fxy(x, y):
                        id = self.xy2id[(x, y)]
                        # 計算特徵函式f(x,y)關於模型P(Y|X)與經驗分佈P~(X)的期望值
                        self.EPx[id] += Pyx * (1.0 / self.N)

    def fxy(self, x, y):
        return (x, y) in self.xy2id

    def train(self, X, Y):
        # 初始化引數，計算一些概率，期望等變數
        self.init_params(X, Y)
        # 初始化每個樣本特徵函式的權重w = 0
        self.w = [0.0 for i in range(self.n)]

        max_iteration = 50
        for times in xrange(max_iteration):
            print 'iterater times %d' % times
            sigmas = []
            self.cal_EPx()

            for i in xrange(self.n):
                # 書上91頁，公式6.34
                sigma = 1 / self.M * math.log(self.EPxy[i] / self.EPx[i])
                sigmas.append(sigma)

            # if len(filter(lambda x: abs(x) >= 0.01, sigmas)) == 0:
            #     break
            # 一次訓練一個樣本，一次更新所有特徵函式的權重w值
            self.w = [self.w[i] + sigmas[i] for i in xrange(self.n)]

    def predict(self, testset):
        results = []
        for test in testset:
            # 基於訓練得到的特徵函式，預測新的待檢測樣本在最大熵模型下的P(Y|x)的argmax(Y)值
            result = self.cal_probality(test)
            results.append(max(result, key=lambda x: x[0])[1])
        return results


def rebuild_features(features):
    '''
    將原feature的（a0,a1,a2,a3,a4,...）
    變成 (0_a0,1_a1,2_a2,3_a3,4_a4,...)形式
    '''
    new_features = []
    for feature in features:
        new_feature = []
        for i, f in enumerate(feature):
            new_feature.append(str(i) + '_' + str(f))
        new_features.append(new_feature)
    return new_features


if __name__ == "__main__":

    print 'Start read data'

    time_1 = time.time()

    path = './data/train_binary.csv'
    pwd = os.getcwd()
    os.chdir(os.path.dirname(path))
    data = pd.read_csv(os.path.basename(path), encoding='gbk')
    os.chdir(pwd)

    raw_data = pd.read_csv(path, header=0)
    data = raw_data.values

    imgs = data[0::, 1::]
    labels = data[::, 0]

    # 選取 2/3 資料作為訓練集， 1/3 資料作為測試集
    train_features, test_features, train_labels, test_labels = train_test_split(
        imgs, labels, test_size=0.33, random_state=23323)

    # 將每個畫素點作為一個特徵維度，構造特徵空間
    train_features = rebuild_features(train_features)
    test_features = rebuild_features(test_features)

    time_2 = time.time()
    print 'read data cost ', time_2 - time_1, ' second', '\n'

    print 'Start training'
    met = MaxEnt()
    met.train(train_features, train_labels)

    time_3 = time.time()
    print 'training cost ', time_3 - time_2, ' second', '\n'

    print 'Start predicting'
    test_predict = met.predict(test_features)
    time_4 = time.time()
    print 'predicting cost ', time_4 - time_3, ' second', '\n'

    score = accuracy_score(test_labels, test_predict)
    print "The accruacy socre is ", score

 

7. 基於scikit-learn實驗logistic regression

0x1：MNIST classfification using multinomial logistic + L1

# -*- coding:utf-8 -*-

import time
import matplotlib.pyplot as plt
import numpy as np

from sklearn.datasets import fetch_mldata
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.utils import check_random_state

# Turn down for faster convergence
t0 = time.time()
train_samples = 5000

mnist = fetch_mldata('MNIST original')
X = mnist.data.astype('float64')
y = mnist.target
random_state = check_random_state(0)
permutation = random_state.permutation(X.shape[0])
X = X[permutation]
y = y[permutation]
X = X.reshape((X.shape[0], -1))

X_train, X_test, y_train, y_test = train_test_split(
    X, y, train_size=train_samples, test_size=10000)

# 特徵標準化
scaler = StandardScaler()
X_train = scaler.fit_transform(X_train)
X_test = scaler.transform(X_test)

# Turn up tolerance for faster convergence
clf = LogisticRegression(C=50. / train_samples,
                         multi_class='multinomial',
                         penalty='l1', solver='saga', tol=0.1)
clf.fit(X_train, y_train)
sparsity = np.mean(clf.coef_ == 0) * 100
score = clf.score(X_test, y_test)
# print('Best C % .4f' % clf.C_)
print("Sparsity with L1 penalty: %.2f%%" % sparsity)
print("Test score with L1 penalty: %.4f" % score)

coef = clf.coef_.copy()
plt.figure(figsize=(10, 5))
scale = np.abs(coef).max()
for i in range(10):
    l1_plot = plt.subplot(2, 5, i + 1)
    l1_plot.imshow(coef[i].reshape(28, 28), interpolation='nearest',
                   cmap=plt.cm.RdBu, vmin=-scale, vmax=scale)
    l1_plot.set_xticks(())
    l1_plot.set_yticks(())
    l1_plot.set_xlabel('Class %i' % i)
plt.suptitle('Classification vector for...')

run_time = time.time() - t0
print('Example run in %.3f s' % run_time)
plt.show()

0x2：展示了邏輯斯蒂迴歸和簡單線性迴歸對非非線性資料集的分類能力

# -*- coding:utf-8 -*-

import numpy as np
import matplotlib.pyplot as plt

from sklearn import linear_model

# this is our test set, it's just a straight line with some
# Gaussian noise
xmin, xmax = -5, 5
n_samples = 100
np.random.seed(0)
X = np.random.normal(size=n_samples)
y = (X > 0).astype(np.float)
X[X > 0] *= 4
X += .3 * np.random.normal(size=n_samples)

X = X[:, np.newaxis]
# run the classifier
clf = linear_model.LogisticRegression(C=1e5)
clf.fit(X, y)

# and plot the result
plt.figure(1, figsize=(4, 3))
plt.clf()
plt.scatter(X.ravel(), y, color='black', zorder=20)
X_test = np.linspace(-5, 10, 300)


def model(x):
    return 1 / (1 + np.exp(-x))
loss = model(X_test * clf.coef_ + clf.intercept_).ravel()
plt.plot(X_test, loss, color='red', linewidth=3)

ols = linear_model.LinearRegression()
ols.fit(X, y)
plt.plot(X_test, ols.coef_ * X_test + ols.intercept_, linewidth=1)
plt.axhline(.5, color='.5')

plt.ylabel('y')
plt.xlabel('X')
plt.xticks(range(-5, 10))
plt.yticks([0, 0.5, 1])
plt.ylim(-.25, 1.25)
plt.xlim(-4, 10)
plt.legend(('Logistic Regression Model', 'Linear Regression Model'),
           loc="lower right", fontsize='small')
plt.show()

Relevant Link:

http://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LogisticRegression.html
http://scikit-learn.org/stable/auto_examples/linear_model/plot_sparse_logistic_regression_mnist.html#sphx-glr-auto-examples-linear-model-plot-sparse-logistic-regression-mnist-py
http://scikit-learn.org/stable/auto_examples/linear_model/plot_logistic.html#sphx-glr-auto-examples-linear-model-plot-logistic-py

Copyright (c) 2017 LittleHann All rights reserved