第6章 逻辑斯蒂回归和最大熵模型

逻辑斯谛回归(LR)是经典的分类方法

1.逻辑斯谛回归模型是由以下条件概率分布表示的分类模型。逻辑斯谛回归模型可以用于二类或多类分类。

P ( Y = k ∣ x ) = exp ⁡ ( w k ⋅ x ) 1 + ∑ k = 1 K − 1 exp ⁡ ( w k ⋅ x ) , k = 1 , 2 , ⋯   , K − 1 P(Y=k | x)=\frac{\exp \left(w_{k} \cdot x\right)}{1+\sum_{k=1}^{K-1} \exp \left(w_{k} \cdot x\right)}, \quad k=1,2, \cdots, K-1 P(Y=kx)=1+k=1K1exp(wkx)exp(wkx),k=1,2,,K1

P ( Y = K ∣ x ) = 1 1 + ∑ k = 1 K − 1 exp ⁡ ( w k ⋅ x ) P(Y=K | x)=\frac{1}{1+\sum_{k=1}^{K-1} \exp \left(w_{k} \cdot x\right)} P(Y=Kx)=1+k=1K1exp(wkx)1

这里, x x x为输入特征, w w w为特征的权值。

逻辑斯谛回归模型源自逻辑斯谛分布,其分布函数 F ( x ) F(x) F(x) S S S形函数。逻辑斯谛回归模型是由输入的线性函数表示的输出的对数几率模型。

2.最大熵模型是由以下条件概率分布表示的分类模型。最大熵模型也可以用于二类或多类分类。
P w ( y ∣ x ) = 1 Z w ( x ) exp ⁡ ( ∑ i = 1 n w i f i ( x , y ) ) P_{w}(y | x)=\frac{1}{Z_{w}(x)} \exp \left(\sum_{i=1}^{n} w_{i} f_{i}(x, y)\right) Pw(yx)=Zw(x)1exp(i=1nwifi(x,y))

Z w ( x ) = ∑ y exp ⁡ ( ∑ i = 1 n w i f i ( x , y ) ) Z_{w}(x)=\sum_{y} \exp \left(\sum_{i=1}^{n} w_{i} f_{i}(x, y)\right) Zw(x)=yexp(i=1nwifi(x,y))

其中, Z w ( x ) Z_w(x) Zw(x)是规范化因子, f i f_i fi为特征函数, w i w_i wi为特征的权值。

3.最大熵模型可以由最大熵原理推导得出。最大熵原理是概率模型学习或估计的一个准则。最大熵原理认为在所有可能的概率模型(分布)的集合中,熵最大的模型是最好的模型。

最大熵原理应用到分类模型的学习中,有以下约束最优化问题:

min ⁡ − H ( P ) = ∑ x , y P ~ ( x ) P ( y ∣ x ) log ⁡ P ( y ∣ x ) \min -H(P)=\sum_{x, y} \tilde{P}(x) P(y | x) \log P(y | x) minH(P)=x,yP~(x)P(yx)logP(yx)

s . t . P ( f i ) − P ~ ( f i ) = 0 , i = 1 , 2 , ⋯   , n s.t. \quad P\left(f_{i}\right)-\tilde{P}\left(f_{i}\right)=0, \quad i=1,2, \cdots, n s.t.P(fi)P~(fi)=0,i=1,2,,n

∑ y P ( y ∣ x ) = 1 \sum_{y} P(y | x)=1 yP(yx)=1

求解此最优化问题的对偶问题得到最大熵模型。

4.逻辑斯谛回归模型与最大熵模型都属于对数线性模型。

5.逻辑斯谛回归模型及最大熵模型学习一般采用极大似然估计,或正则化的极大似然估计。逻辑斯谛回归模型及最大熵模型学习可以形式化为无约束最优化问题。求解该最优化问题的算法有改进的迭代尺度法、梯度下降法、拟牛顿法。

这里也感谢DataWhale的无私奉献,我结合了原先的学习资料,整合了一下,得到现有的笔记,这里也欢迎大家参与DataWhale的组队学习中,https://datawhalechina.github.io/statistical-learning-method-solutions-manual/#/


最大熵模型

  • 最大熵模型(maximun entropy model)由最大熵原理推导实现,而最大熵原理损失概率模型学习的一个准则。最大熵原理认为,学习概率模型是,所有可能的概率模型(分布)中,熵最大的模型是最好的模型。

  • 最大熵原理也可以表述为在满足约束条件的模型集合中取熵最大的模型

  • 假设满足所有约束条件的模型的集合为
    KaTeX parse error: Got function '\tilde' with no arguments as subscript at position 29: …n P|E_P(f_i)=E_\̲t̲i̲l̲d̲e̲{P}({f_i})\}

  • 定义在条件概率分布 P ( X ∣ Y ) P(X|Y) P(XY)上的条件熵为
    H ( P ) = − ∑ x , y P ~ ( x ) P ( y ∣ x ) log ⁡ P ( y ∣ x ) H(P)=-\sum_{x,y}\tilde{P}(x)P(y|x)\log P(y|x) H(P)=x,yP~(x)P(yx)logP(yx)
    则模型集合C中条件熵 H ( P ) H(P) H(P)最大的模型称为最大熵模型。式中的对数为自然对数。

最大熵模型的学习

  • 最大熵模型的学习过程就是求解最大熵模型的过程。最大熵模型的学习可以形式化为约束最优的问题。最大熵模型的学习就等价与约束最优化问题
  • 将约束最优化的原始问题转换为无约束最优化的对偶问题。通过求解对偶问题求解原始问题。
  • 简单来说,约束最优化问题包含 ≤ 0 \leq 0 0,和 = 0 =0 =0两种约束条件,这是约束优化问题的一般形式

min ⁡ x ∈ R n f ( x ) s . t . c i ( x ) ≤ 0 , i = 1 , 2 , … , k h j ( x ) = 0 , j = 1 , 2 , … , l \begin{aligned} \min_{x \in R^n}\quad &f(x) \\ s.t.\quad&c_i(x) \leq 0 , i=1,2,\ldots,k\\ &h_j(x) = 0 , j=1,2,\ldots,l \end{aligned} xRnmins.t.f(x)ci(x)0,i=1,2,,khj(x)=0,j=1,2,,l

  • 引入广义拉格朗日函数

L ( x , α , β ) = f ( x ) + ∑ i = 0 k α i c i ( x ) + ∑ j = 1 l β j h j ( x ) L(x,\alpha,\beta) = f(x) + \sum_{i=0}^k \alpha_ic_i(x) + \sum_{j=1}^l \beta_jh_j(x) L(x,α,β)=f(x)+i=0kαici(x)+j=1lβjhj(x)

α \alpha α β \beta β是拉格朗日乘子

  • 在KKT的条件下,原始问题和对偶问题的最优值相等

∇ x L ( x ∗ , α ∗ , β ∗ ) = 0 ∇ α L ( x ∗ , α ∗ , β ∗ ) = 0 ∇ β L ( x ∗ , α ∗ , β ∗ ) = 0 α i ∗ c i ( x ∗ ) = 0 , i = 1 , 2 , … , k c i ( x ∗ ) ≤ 0 , i = 1 , 2 , … , k α i ∗ ≥ 0 , i = 1 , 2 , … , k h j ( x ∗ ) = 0 , j = 1 , 2 , … , l ∇_xL(x^∗,α^∗,β^∗)=0\\ ∇_αL(x^∗,α^∗,β^∗)=0\\ ∇_βL(x^∗,α^∗,β^∗)=0\\ α_i^∗c_i(x^*)=0,i=1,2,…,k\\ c_i(x^*)≤0,i=1,2,…,k\\ α^∗_i≥0,i=1,2,…,k\\ h_j(x^∗)=0,j=1,2,…,l xL(x,α,β)=0αL(x,α,β)=0βL(x,α,β)=0αici(x)=0,i=1,2,,kci(x)0,i=1,2,,kαi0,i=1,2,,khj(x)=0,j=1,2,,l

  • 前面三个条件是由解析函数的知识,对于各个变量的偏导数为0,后面四个条件就是原始问题的约束条件以及拉格朗日乘子需要满足的约束,第四个条件是KKT的对偶互补条件

  • L ( P , w ) L(P, w) L(P,w) P P P求偏导并令其为零解得

P ( y ∣ x ) = exp ⁡ ( ∑ i = 1 n w i f i ( x , y ) + w 0 − 1 ) = exp ⁡ ( ∑ i = 1 n w i f i ( x , y ) ) exp ⁡ ( 1 − w 0 ) P(y|x)=\exp{\left(\sum_{i=1}^{n}w_if_i(x,y)+w_0-1\right)}=\frac{\exp{\left(\sum\limits_{i=1}^{n}w_if_i(x,y)\right)}}{\exp{\left(1-w_0\right)}} P(yx)=exp(i=1nwifi(x,y)+w01)=exp(1w0)exp(i=1nwifi(x,y))

  • 因为 ∑ y P ( y ∣ x ) = 1 \sum\limits_{y}P(y|x)=1 yP(yx)=1,然后得到模型

P w ( y ∣ x ) = 1 Z w ( x ) exp ⁡ ∑ i = 1 n w i f i ( x , y ) P_w(y|x)=\frac{1}{Z_w(x)}\exp{\sum\limits_{i=1}^{n}w_if_i(x,y)}\\ Pw(yx)=Zw(x)1expi=1nwifi(x,y)

其中, Z w ( x ) = ∑ y exp ⁡ ( ∑ i = 1 n w i f i ( x , y ) ) 其中,Z_w(x)=\sum_{y}\exp({\sum_{i=1}^{n}w_if_i(x,y))} 其中,Zw(x)=yexp(i=1nwifi(x,y))

  • 这里 Z w ( x ) Z_w(x) Zw(x)先用来代替 exp ⁡ ( 1 − w 0 ) \exp(1-w_0) exp(1w0), Z w Z_w Zw是归一化因子
  • 并不是因为概率为1推导出了 Z w Z_w Zw的表达式,这样一个表达式是凑出来的,意思就是遍历 y y y的所有取值,求分子表达式的占比

极大似然估计

  • 对偶函数的极大化等价于最大熵模型的极大似然估计
  • 已知训练数据的经验分布 P ~ ( X , Y ) \widetilde {P}(X,Y) P (X,Y),条件概率分布 P ( Y ∣ X ) P(Y|X) P(YX)的对数似然函数表示为

L P ~ ( P w ) = log ⁡ ∏ x , y P ( y ∣ x ) P ~ ( x , y ) = ∑ x , y P ~ ( x , y ) log ⁡ P ( y ∣ x ) L_{\widetilde {P}}(P_w)=\log\prod_{x,y}P(y|x)^{\widetilde {P}(x,y)}=\sum \limits_{x,y}\widetilde {P}(x,y)\log{P}(y|x) LP (Pw)=logx,yP(yx)P (x,y)=x,yP (x,y)logP(yx)

  • 当条件分布概率 P ( y ∣ x ) P(y|x) P(yx)是最大熵模型时

L P ~ ( P w ) = ∑ x , y P ~ ( x , y ) log ⁡ P ( y ∣ x ) = ∑ x , y P ~ ( x , y ) ∑ i = 1 n w i f i ( x , y ) − ∑ x , y P ~ ( x , y ) log ⁡ ( Z w ( x ) ) = ∑ x , y P ~ ( x , y ) ∑ i = 1 n w i f i ( x , y ) − ∑ x , y P ~ ( x ) P ( y ∣ x ) log ⁡ ( Z w ( x ) ) = ∑ x , y P ~ ( x , y ) ∑ i = 1 n w i f i ( x , y ) − ∑ x P ~ ( x ) log ⁡ ( Z w ( x ) ) ∑ y P ( y ∣ x ) = ∑ x , y P ~ ( x , y ) ∑ i = 1 n w i f i ( x , y ) − ∑ x P ~ ( x ) log ⁡ ( Z w ( x ) )   \begin{aligned} L_{\widetilde {P}}(P_w)&=\sum \limits_{x,y}\widetilde {P}(x,y)\log{P}(y|x)\\ &=\sum \limits_{x,y}\widetilde {P}(x,y)\sum \limits_{i=1}^{n}w_if_i(x,y) -\sum \limits_{x,y}\widetilde{P}(x,y)\log{(Z_w(x))}\\ &=\sum \limits_{x,y}\widetilde {P}(x,y)\sum \limits_{i=1}^{n}w_if_i(x,y) -\sum \limits_{x,y}\widetilde{P}(x)P(y|x)\log{(Z_w(x))}\\ &=\sum \limits_{x,y}\widetilde {P}(x,y)\sum \limits_{i=1}^{n}w_if_i(x,y) -\sum \limits_{x}\widetilde{P}(x)\log{(Z_w(x))}\sum_{y}P(y|x)\\ &=\sum \limits_{x,y}\widetilde {P}(x,y)\sum \limits_{i=1}^{n}w_if_i(x,y) -\sum \limits_{x}\widetilde{P}(x)\log{(Z_w(x))} \end{aligned} \ LP (Pw)=x,yP (x,y)logP(yx)=x,yP (x,y)i=1nwifi(x,y)x,yP (x,y)log(Zw(x))=x,yP (x,y)i=1nwifi(x,y)x,yP (x)P(yx)log(Zw(x))=x,yP (x,y)i=1nwifi(x,y)xP (x)log(Zw(x))yP(yx)=x,yP (x,y)i=1nwifi(x,y)xP (x)log(Zw(x)) 

  • 推导过程用到了 ∑ y P ( y ∣ x ) = 1 \sum\limits_yP(y|x)=1 yP(yx)=1

模型学习的最优化算法

  • 逻辑斯蒂回归模型、最大熵模型学习归结为以似然函数为目标函数的最优化问题,通常通过迭代算法求解
  • 用最优化的观点来看,常用找最优解的方法有改建的迭代尺度法、梯度下降法、牛顿法或拟牛顿法。牛顿法和拟牛顿法一般收敛速度更快。
  • 简单算法在这里不做描述了,因为之前也做了笔记,这里着重描述一下改进的迭代尺度算法

改进的迭代尺度算法(IIS)

输入:特征函数 f 1 , f 2 , … , f n f_1,f_2,\dots,f_n f1,f2,,fn;经验分布 P ~ ( X , Y ) \tilde{P}(X,Y) P~(X,Y),模型 P w ( y ∣ x ) P_w(y|x) Pw(yx);

输出:最优参数值 w i ∗ w_i^* wi;最优模型 P w ∗ P_{w^*} Pw

(1)对所有 i ∈ { 1 , 2 , … , n } i \in \{ 1,2,\dots,n\} i{1,2,,n},取初值 w i = 0 w_i=0 wi=0

(2)对每一 i ∈ { 1 , 2 , … , n } i \in \{ 1,2,\dots,n\} i{1,2,,n}

(a) 令 δ i \delta_{i} δi 是方程
∑ x , y P ~ ( x ) P ( y ∣ x ) f i ( x , y ) exp ⁡ ( δ i f # ( x , y ) ) = E P ~ ( f i ) \sum_{x, y} \tilde{P}(x) P(y \mid x) f_{i}(x, y) \exp \left(\delta_{i} f^{\#}(x, y)\right)=E_{\tilde{P}}\left(f_{i}\right) x,yP~(x)P(yx)fi(x,y)exp(δif#(x,y))=EP~(fi)
的解, 这里,
f # ( x , y ) = ∑ i = 1 n f i ( x , y ) f^{\#}(x, y)=\sum_{i=1}^{n} f_{i}(x, y) f#(x,y)=i=1nfi(x,y)
(b)更新 w i w_{i} wi 值: w i ← w i + δ i w_{i} \leftarrow w_{i}+\delta_{i} wiwi+δi
(3)如果不是所有 w i w_{i} wi 都收玫,重复步 (2)。

这一算法关键的一步是 ( a ) (\mathrm{a}) (a), 即求解方程中的 δ i \delta_{i} δi 。如果 f # ( x , y ) f^{\#}(x, y) f#(x,y) 是常数, 即 对任何 x , y x, y x,y, 有 f # ( x , y ) = M f^{\#}(x, y)=M f#(x,y)=M, 那么 δ i \delta_{i} δi 可以显式地表示成
δ i = 1 M log ⁡ E P ~ ( f i ) E P ( f i ) \delta_{i}=\frac{1}{M} \log \frac{E_{\tilde{P}}\left(f_{i}\right)}{E_{P}\left(f_{i}\right)} δi=M1logEP(fi)EP~(fi)
如果 f # ( x , y ) f^{\#}(x, y) f#(x,y) 不是常数,那么必须通过数值计算求 δ i \delta_{i} δi 。简单有效的方法是牛顿法。

g ( δ i ) = 0 g\left(\delta_{i}\right)=0 g(δi)=0 表示方程, 牛顿法通过迭代求得 δ i ∗ \delta_{i}^{*} δi, 使得 g ( δ i ∗ ) = 0 g\left(\delta_{i}^{*}\right)=0 g(δi)=0 。迭代公式是
δ i ( k + 1 ) = δ i ( k ) − g ( δ i ( k ) ) g ′ ( δ i ( k ) ) \delta_{i}^{(k+1)}=\delta_{i}^{(k)}-\frac{g\left(\delta_{i}^{(k)}\right)}{g^{\prime}\left(\delta_{i}^{(k)}\right)} δi(k+1)=δi(k)g(δi(k))g(δi(k))
只要适当选取初始值 δ i ( 0 ) \delta_{i}^{(0)} δi(0), 由于 δ i \delta_{i} δi 的方程有单根, 因此牛顿法恒收敛, 而且收敛速度很快。

练习

class MaxEntropy:
def init(self, EPS=0.005):
self._samples = []
self._Y = set() # 标签集合,相当去去重后的y
self._numXY = {} # key为(x,y),value为出现次数
self._N = 0 # 样本数
self.Ep = [] # 样本分布的特征期望值
self._xyID = {} # key记录(x,y),value记录id号
self._n = 0 # 特征键值(x,y)的个数
self._C = 0 # 最大特征数
self._IDxy = {} # key为(x,y),value为对应的id号
self._w = []
self._EPS = EPS # 收敛条件
self._lastw = [] # 上一次w参数值

def loadData(self, dataset):
    self._samples = deepcopy(dataset)
    for items in self._samples:
        y = items[0]
        X = items[1:]
        self._Y.add(y)  # 集合中y若已存在则会自动忽略
        for x in X:
            if (x, y) in self._numXY:
                self._numXY[(x, y)] += 1
            else:
                self._numXY[(x, y)] = 1

    self._N = len(self._samples)
    self._n = len(self._numXY)
    self._C = max([len(sample) - 1 for sample in self._samples])
    self._w = [0] * self._n
    self._lastw = self._w[:]

    self._Ep_ = [0] * self._n
    for i, xy in enumerate(self._numXY):  # 计算特征函数fi关于经验分布的期望
        self._Ep_[i] = self._numXY[xy] / self._N
        self._xyID[xy] = i
        self._IDxy[i] = xy

def _Zx(self, X):  # 计算每个Z(x)值
    zx = 0
    for y in self._Y:
        ss = 0
        for x in X:
            if (x, y) in self._numXY:
                ss += self._w[self._xyID[(x, y)]]
        zx += math.exp(ss)
    return zx

def _model_pyx(self, y, X):  # 计算每个P(y|x)
    zx = self._Zx(X)
    ss = 0
    for x in X:
        if (x, y) in self._numXY:
            ss += self._w[self._xyID[(x, y)]]
    pyx = math.exp(ss) / zx
    return pyx

def _model_ep(self, index):  # 计算特征函数fi关于模型的期望
    x, y = self._IDxy[index]
    ep = 0
    for sample in self._samples:
        if x not in sample:
            continue
        pyx = self._model_pyx(y, sample)
        ep += pyx / self._N
    return ep

def _convergence(self):  # 判断是否全部收敛
    for last, now in zip(self._lastw, self._w):
        if abs(last - now) >= self._EPS:
            return False
    return True

def predict(self, X):  # 计算预测概率
    Z = self._Zx(X)
    result = {}
    for y in self._Y:
        ss = 0
        for x in X:
            if (x, y) in self._numXY:
                ss += self._w[self._xyID[(x, y)]]
        pyx = math.exp(ss) / Z
        result[y] = pyx
    return result

def train(self, maxiter=1000):  # 训练数据
    for loop in range(maxiter):  # 最大训练次数
        print("iter:%d" % loop)
        self._lastw = self._w[:]
        for i in range(self._n):
            ep = self._model_ep(i)  # 计算第i个特征的模型期望
            self._w[i] += math.log(self._Ep_[i] / ep) / self._C  # 更新参数
        print("w:", self._w)
        if self._convergence():  # 判断是否收敛
            break
dataset = [['no', 'sunny', 'hot', 'high', 'FALSE'],
           ['no', 'sunny', 'hot', 'high', 'TRUE'],
           ['yes', 'overcast', 'hot', 'high', 'FALSE'],
           ['yes', 'rainy', 'mild', 'high', 'FALSE'],
           ['yes', 'rainy', 'cool', 'normal', 'FALSE'],
           ['no', 'rainy', 'cool', 'normal', 'TRUE'],
           ['yes', 'overcast', 'cool', 'normal', 'TRUE'],
           ['no', 'sunny', 'mild', 'high', 'FALSE'],
           ['yes', 'sunny', 'cool', 'normal', 'FALSE'],
           ['yes', 'rainy', 'mild', 'normal', 'FALSE'],
           ['yes', 'sunny', 'mild', 'normal', 'TRUE'],
           ['yes', 'overcast', 'mild', 'high', 'TRUE'],
           ['yes', 'overcast', 'hot', 'normal', 'FALSE'],
           ['no', 'rainy', 'mild', 'high', 'TRUE']]
maxent = MaxEntropy()
x = ['overcast', 'mild', 'high', 'FALSE']
maxent.loadData(dataset)
maxent.train()
print('predict:', maxent.predict(x))
predict: {'yes': 0.9999971802186581, 'no': 2.819781341881656e-06}

第6章Logistic回归与最大熵模型-习题

习题6.1

​ 确认Logistic分布属于指数分布族。

解答:
第1步:
首先给出指数分布族的定义:
对于随机变量 x x x,在给定参数 η \eta η下,其概率分别满足如下形式: p ( x ∣ η ) = h ( x ) g ( η ) exp ⁡ ( η T u ( x ) ) p(x|\eta)=h(x)g(\eta)\exp(\eta^Tu(x)) p(xη)=h(x)g(η)exp(ηTu(x))我们称之为指数分布族
其中:
x x x:可以是标量或者向量,可以是离散值也可以是连续值
η \eta η:自然参数
g ( η ) g(\eta) g(η):归一化系数
h ( x ) , u ( x ) h(x),u(x) h(x),u(x) x x x的某个函数


**第2步:**证明伯努利分布属于指数分布族
伯努利分布: φ \varphi φ y = 1 y=1 y=1的概率,即 P ( Y = 1 ) = φ P(Y=1)=\varphi P(Y=1)=φ
P ( y ∣ φ ) = φ y ( 1 − φ ) ( 1 − y ) = ( 1 − φ ) φ y ( 1 − φ ) ( − y ) = ( 1 − φ ) ( φ 1 − φ ) y = ( 1 − φ ) exp ⁡ ( y ln ⁡ φ 1 − φ ) = 1 1 + e η exp ⁡ ( η y ) \begin{aligned} P(y|\varphi) &= \varphi^y (1-\varphi)^{(1-y)} \\ &= (1-\varphi) \varphi^y (1-\varphi)^{(-y)} \\ &= (1-\varphi) (\frac{\varphi}{1-\varphi})^y \\ &= (1-\varphi) \exp\left(y \ln \frac{\varphi}{1-\varphi} \right) \\ &= \frac{1}{1+e^\eta} \exp (\eta y) \end{aligned} P(yφ)=φy(1φ)(1y)=(1φ)φy(1φ)(y)=(1φ)(1φφ)y=(1φ)exp(yln1φφ)=1+eη1exp(ηy)
其中, η = ln ⁡ φ 1 − φ ⇔ φ = 1 1 + e − η \displaystyle \eta=\ln \frac{\varphi}{1-\varphi} \Leftrightarrow \varphi = \frac{1}{1+e^{-\eta}} η=ln1φφφ=1+eη1
y y y替换成 x x x,可得 P ( x ∣ η ) = 1 1 + e η exp ⁡ ( η x ) \displaystyle P(x|\eta) = \frac{1}{1+e^\eta} \exp (\eta x) P(xη)=1+eη1exp(ηx)
对比可知,伯努利分布属于指数分布族,其中 h ( x ) = 1 , g ( η ) = 1 1 + e η , u ( x ) = x \displaystyle h(x) = 1, g(\eta)= \frac{1}{1+e^\eta}, u(x)=x h(x)=1,g(η)=1+eη1,u(x)=x


第3步:
广义线性模型(GLM)必须满足三个假设:

  1. y ∣ x ; θ ∼ E x p o n e n t i a l F a m i l y ( η ) y | x;\theta \sim ExponentialFamily(\eta) yx;θExponentialFamily(η),即假设预测变量 y y y在给定 x x x,以 θ \theta θ为参数的条件概率下,属于以 η \eta η作为自然参数的指数分布族;
  2. 给定 x x x,求解出以 x x x为条件的 T ( y ) T(y) T(y)的期望 E [ T ( y ) ∣ x ] E[T(y)|x] E[T(y)x],即算法输出为 h ( x ) = E [ T ( y ) ∣ x ] h(x)=E[T(y)|x] h(x)=E[T(y)x]
  3. 满足 η = θ T x \eta=\theta^T x η=θTx,即自然参数和输入特征向量 x x x之间线性相关,关系由 θ \theta θ决定,仅当 η \eta η是实数时才有意义,若 η \eta η是一个向量,则 η i = θ i T x \eta_i=\theta_i^T x ηi=θiTx

**第4步:**推导伯努利分布的GLM
已知伯努利分布属于指数分布族,对给定的 x , η x,\eta x,η,求解期望: h θ ( x ) = E [ y ∣ x ; θ ] = 1 ⋅ p ( y = 1 ) + 0 ⋅ p ( y = 0 ) = φ = 1 1 + e − η = 1 1 + e − θ T x \begin{aligned} h_{\theta}(x) &= E[y|x;\theta] \\ &= 1 \cdot p(y=1)+ 0 \cdot p(y=0) \\ &= \varphi \\ &= \frac{1}{1+e^{-\eta}} \\ &= \frac{1}{1+e^{-\theta^T x}} \end{aligned} hθ(x)=E[yx;θ]=1p(y=1)+0p(y=0)=φ=1+eη1=1+eθTx1可得到Logistic回归算法,故Logistic分布属于指数分布族,得证。

习题6.2

  写出Logistic回归模型学习的梯度下降算法。

解答:
对于Logistic模型: P ( Y = 1 ∣ x ) = exp ⁡ ( w ⋅ x + b ) 1 + exp ⁡ ( w ⋅ x + b ) P ( Y = 0 ∣ x ) = 1 1 + exp ⁡ ( w ⋅ x + b ) P(Y=1 | x)=\frac{\exp (w \cdot x+b)}{1+\exp (w \cdot x+b)} \\ P(Y=0 | x)=\frac{1}{1+\exp (w \cdot x+b)} P(Y=1∣x)=1+exp(wx+b)exp(wx+b)P(Y=0∣x)=1+exp(wx+b)1对数似然函数为: L ( w ) = ∑ i = 1 N [ y i ( w ⋅ x i ) − log ⁡ ( 1 + exp ⁡ ( w ⋅ x i ) ) ] \displaystyle L(w)=\sum_{i=1}^N \left[y_i (w \cdot x_i)-\log \left(1+\exp (w \cdot x_i)\right)\right] L(w)=i=1N[yi(wxi)log(1+exp(wxi))]
似然函数求偏导,可得 ∂ L ( w ) ∂ w ( j ) = ∑ i = 1 N [ x i ( j ) ⋅ y i − exp ⁡ ( w ⋅ x i ) ⋅ x i ( j ) 1 + exp ⁡ ( w ⋅ x i ) ] \displaystyle \frac{\partial L(w)}{\partial w^{(j)}}=\sum_{i=1}^N\left[x_i^{(j)} \cdot y_i-\frac{\exp (w \cdot x_i) \cdot x_i^{(j)}}{1+\exp (w \cdot x_i)}\right] w(j)L(w)=i=1N[xi(j)yi1+exp(wxi)exp(wxi)xi(j)]
梯度函数为: ∇ L ( w ) = [ ∂ L ( w ) ∂ w ( 0 ) , ⋯   , ∂ L ( w ) ∂ w ( m ) ] \displaystyle \nabla L(w)=\left[\frac{\partial L(w)}{\partial w^{(0)}}, \cdots, \frac{\partial L(w)}{\partial w^{(m)}}\right] L(w)=[w(0)L(w),,w(m)L(w)]
Logistic回归模型学习的梯度下降算法:
(1) 取初始值 x ( 0 ) ∈ R x^{(0)} \in R x(0)R,置 k = 0 k=0 k=0
(2) 计算 f ( x ( k ) ) f(x^{(k)}) f(x(k))
(3) 计算梯度 g k = g ( x ( k ) ) g_k=g(x^{(k)}) gk=g(x(k)),当 ∥ g k ∥ < ε \|g_k\| < \varepsilon gk<ε时,停止迭代,令 x ∗ = x ( k ) x^* = x^{(k)} x=x(k);否则,求 λ k \lambda_k λk,使得 f ( x ( k ) + λ k g k ) = max ⁡ λ ⩾ 0 f ( x ( k ) + λ g k ) \displaystyle f(x^{(k)}+\lambda_k g_k) = \max_{\lambda \geqslant 0}f(x^{(k)}+\lambda g_k) f(x(k)+λkgk)=λ0maxf(x(k)+λgk)
(4) 置 x ( k + 1 ) = x ( k ) + λ k g k x^{(k+1)}=x^{(k)}+\lambda_k g_k x(k+1)=x(k)+λkgk,计算 f ( x ( k + 1 ) ) f(x^{(k+1)}) f(x(k+1)),当 ∥ f ( x ( k + 1 ) ) − f ( x ( k ) ) ∥ < ε \|f(x^{(k+1)}) - f(x^{(k)})\| < \varepsilon f(x(k+1))f(x(k))<ε ∥ x ( k + 1 ) − x ( k ) ∥ < ε \|x^{(k+1)} - x^{(k)}\| < \varepsilon x(k+1)x(k)<ε时,停止迭代,令 x ∗ = x ( k + 1 ) x^* = x^{(k+1)} x=x(k+1)
(5) 否则,置 k = k + 1 k=k+1 k=k+1,转(3)

%matplotlib inline
import numpy as np
import time
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from pylab import mpl

# 图像显示中文
mpl.rcParams['font.sans-serif'] = ['Microsoft YaHei']


class LogisticRegression:
    def __init__(self, learn_rate=0.1, max_iter=10000, tol=1e-2):
        self.learn_rate = learn_rate  # 学习率
        self.max_iter = max_iter  # 迭代次数
        self.tol = tol  # 迭代停止阈值
        self.w = None  # 权重

    def preprocessing(self, X):
        """将原始X末尾加上一列,该列数值全部为1"""
        row = X.shape[0]
        y = np.ones(row).reshape(row, 1)
        X_prepro = np.hstack((X, y))
        return X_prepro

    def sigmod(self, x):
        return 1 / (1 + np.exp(-x))

    def fit(self, X_train, y_train):
        X = self.preprocessing(X_train)
        y = y_train.T
        # 初始化权重w
        self.w = np.array([[0] * X.shape[1]], dtype=np.float)
        k = 0
        for loop in range(self.max_iter):
            # 计算梯度
            z = np.dot(X, self.w.T)
            grad = X * (y - self.sigmod(z))
            grad = grad.sum(axis=0)
            # 利用梯度的绝对值作为迭代中止的条件
            if (np.abs(grad) <= self.tol).all():
                break
            else:
                # 更新权重w 梯度上升——求极大值
                self.w += self.learn_rate * grad
                k += 1
        print("迭代次数:{}次".format(k))
        print("最终梯度:{}".format(grad))
        print("最终权重:{}".format(self.w[0]))

    def predict(self, x):
        p = self.sigmod(np.dot(self.preprocessing(x), self.w.T))
        print("Y=1的概率被估计为:{:.2%}".format(p[0][0]))  # 调用score时,注释掉
        p[np.where(p > 0.5)] = 1
        p[np.where(p < 0.5)] = 0
        return p

    def score(self, X, y):
        y_c = self.predict(X)
        error_rate = np.sum(np.abs(y_c - y.T)) / y_c.shape[0]
        return 1 - error_rate

    def draw(self, X, y):
        # 分离正负实例点
        y = y[0]
        X_po = X[np.where(y == 1)]
        X_ne = X[np.where(y == 0)]
        # 绘制数据集散点图
        ax = plt.axes(projection='3d')
        x_1 = X_po[0, :]
        y_1 = X_po[1, :]
        z_1 = X_po[2, :]
        x_2 = X_ne[0, :]
        y_2 = X_ne[1, :]
        z_2 = X_ne[2, :]
        ax.scatter(x_1, y_1, z_1, c="r", label="正实例")
        ax.scatter(x_2, y_2, z_2, c="b", label="负实例")
        ax.legend(loc='best')
        # 绘制p=0.5的区分平面
        x = np.linspace(-3, 3, 3)
        y = np.linspace(-3, 3, 3)
        x_3, y_3 = np.meshgrid(x, y)
        a, b, c, d = self.w[0]
        z_3 = -(a * x_3 + b * y_3 + d) / c
        ax.plot_surface(x_3, y_3, z_3, alpha=0.5)  # 调节透明度
        plt.show()
# 训练数据集
X_train = np.array([[3, 3, 3], [4, 3, 2], [2, 1, 2], [1, 1, 1], [-1, 0, 1],
                    [2, -2, 1]])
y_train = np.array([[1, 1, 1, 0, 0, 0]])
# 构建实例,进行训练
clf = LogisticRegression()
clf.fit(X_train, y_train)
clf.draw(X_train, y_train)
迭代次数:3232次
最终梯度:[ 0.00144779  0.00046133  0.00490279 -0.00999848]
最终权重:[  2.96908597   1.60115396   5.04477438 -13.43744079]

在这里插入图片描述

习题6.3

  写出最大熵模型学习的DFP算法。(关于一般的DFP算法参见附录B)

解答:
第1步:
最大熵模型为:
max ⁡ H ( p ) = − ∑ x , y P ( x ) P ( y ∣ x ) log ⁡ P ( y ∣ x ) st. E p ( f i ) − E p ^ ( f i ) = 0 , i = 1 , 2 , ⋯   , n ∑ y P ( y ∣ x ) = 1 \begin{array}{cl} {\max } & {H(p)=-\sum_{x, y} P(x) P(y | x) \log P(y | x)} \\ {\text {st.}} & {E_p(f_i)-E_{\hat{p}}(f_i)=0, \quad i=1,2, \cdots,n} \\ & {\sum_y P(y | x)=1} \end{array} maxst.H(p)=x,yP(x)P(yx)logP(yx)Ep(fi)Ep^(fi)=0,i=1,2,,nyP(yx)=1
引入拉格朗日乘子,定义拉格朗日函数:
L ( P , w ) = ∑ x y P ( x ) P ( y ∣ x ) log ⁡ P ( y ∣ x ) + w 0 ( 1 − ∑ y P ( y ∣ x ) ) + ∑ i = 1 w i ( ∑ x y P ( x , y ) f i ( x , y ) − ∑ x y P ( x , y ) P ( y ∣ x ) f i ( x , y ) ) L(P, w)=\sum_{xy} P(x) P(y | x) \log P(y | x)+w_0 \left(1-\sum_y P(y | x)\right) \\ +\sum_{i=1} w_i\left(\sum_{xy} P(x, y) f_i(x, y)-\sum_{xy} P(x, y) P(y | x) f_i(x, y)\right) L(P,w)=xyP(x)P(yx)logP(yx)+w0(1yP(yx))+i=1wi(xyP(x,y)fi(x,y)xyP(x,y)P(yx)fi(x,y))
$
最优化原始问题为:$
min ⁡ P ∈ C max ⁡ w L ( P , w ) \min_{P \in C} \max_{w} L(P,w) PCminwmaxL(P,w)
对偶问题为:
max ⁡ w min ⁡ P ∈ C L ( P , w ) \max_{w} \min_{P \in C} L(P,w) wmaxPCminL(P,w)

Ψ ( w ) = min ⁡ P ∈ C L ( P , w ) = L ( P w , w ) \Psi(w) = \min_{P \in C} L(P,w) = L(P_w, w) Ψ(w)=PCminL(P,w)=L(Pw,w)
Ψ ( w ) \Psi(w) Ψ(w)称为对偶函数,同时,其解记作
P w = arg ⁡ min ⁡ P ∈ C L ( P , w ) = P w ( y ∣ x ) P_w = \mathop{\arg \min}_{P \in C} L(P,w) = P_w(y|x) Pw=argminPCL(P,w)=Pw(yx)
L ( P , w ) L(P,w) L(P,w) P ( y ∣ x ) P(y|x) P(yx)的偏导数,并令偏导数等于0,解得:
P w ( y ∣ x ) = 1 Z w ( x ) exp ⁡ ( ∑ i = 1 n w i f i ( x , y ) ) P_w(y | x)=\frac{1}{Z_w(x)} \exp \left(\sum_{i=1}^n w_i f_i (x, y)\right) Pw(yx)=Zw(x)1exp(i=1nwifi(x,y))
其中:
Z w ( x ) = ∑ y exp ⁡ ( ∑ i = 1 n w i f i ( x , y ) ) Z_w(x)=\sum_y \exp \left(\sum_{i=1}^n w_i f_i(x, y)\right) Zw(x)=yexp(i=1nwifi(x,y))
则最大熵模型目标函数表示为​
φ ( w ) = min ⁡ w ∈ R n Ψ ( w ) = ∑ x P ( x ) log ⁡ ∑ y exp ⁡ ( ∑ i = 1 n w i f i ( x , y ) ) − ∑ x , y P ( x , y ) ∑ i = 1 n w i f i ( x , y ) \varphi(w)=\min_{w \in R_n} \Psi(w) = \sum_{x} P(x) \log \sum_{y} \exp \left(\sum_{i=1}^{n} w_{i} f_{i}(x, y)\right)-\sum_{x, y} P(x, y) \sum_{i=1}^{n} w_{i} f_{i}(x, y) φ(w)=wRnminΨ(w)=xP(x)logyexp(i=1nwifi(x,y))x,yP(x,y)i=1nwifi(x,y)
第2步:
DFP的 G k + 1 G_{k+1} Gk+1的迭代公式为:
G k + 1 = G k + δ k δ k T δ k T y k − G k y k y k T G k y k T G k y k G_{k+1}=G_k+\frac{\delta_k \delta_k^T}{\delta_k^T y_k}-\frac{G_k y_k y_k^T G_k}{y_k^T G_k y_k} Gk+1=Gk+δkTykδkδkTykTGkykGkykykTGk

最大熵模型的DFP算法:
输入:目标函数 φ ( w ) \varphi(w) φ(w),梯度 g ( w ) = ∇ g ( w ) g(w) = \nabla g(w) g(w)=g(w),精度要求 ε \varepsilon ε
输出: φ ( w ) \varphi(w) φ(w)的极小值点 w ∗ w^* w
(1)选定初始点 w ( 0 ) w^{(0)} w(0),取 G 0 G_0 G0为正定对称矩阵,置 k = 0 k=0 k=0
(2)计算 g k = g ( w ( k ) ) g_k=g(w^{(k)}) gk=g(w(k)),若 ∥ g k ∥ < ε \|g_k\| < \varepsilon gk<ε,则停止计算,得近似解 w ∗ = w ( k ) w^*=w^{(k)} w=w(k),否则转(3)
(3)置 p k = − G k g k p_k=-G_kg_k pk=Gkgk
(4)一维搜索:求 λ k \lambda_k λk使得 φ ( w ( k ) + λ k P k ) = min ⁡ λ ⩾ 0 φ ( w ( k ) + λ P k ) \varphi\left(w^{(k)}+\lambda_k P_k\right)=\min _{\lambda \geqslant 0} \varphi\left(w^{(k)}+\lambda P_{k}\right) φ(w(k)+λkPk)=λ0minφ(w(k)+λPk)(5)置 w ( k + 1 ) = w ( k ) + λ k p k w^{(k+1)}=w^{(k)}+\lambda_k p_k w(k+1)=w(k)+λkpk
(6)计算 g k + 1 = g ( w ( k + 1 ) ) g_{k+1}=g(w^{(k+1)}) gk+1=g(w(k+1)),若 ∥ g k + 1 ∥ < ε \|g_{k+1}\| < \varepsilon gk+1<ε,则停止计算,得近似解 w ∗ = w ( k + 1 ) w^*=w^{(k+1)} w=w(k+1);否则,按照迭代式算出 G k + 1 G_{k+1} Gk+1
(7)置 k = k + 1 k=k+1 k=k+1,转(3)

Logo

魔乐社区(Modelers.cn) 是一个中立、公益的人工智能社区,提供人工智能工具、模型、数据的托管、展示与应用协同服务,为人工智能开发及爱好者搭建开放的学习交流平台。社区通过理事会方式运作,由全产业链共同建设、共同运营、共同享有,推动国产AI生态繁荣发展。

更多推荐