最大熵模型

内容源于公众号《简博士数据分析吧》和b站《【合集】十分钟 机器学习 系列视频 《统计学习方法》》
代码来自于https://blog.csdn.net/weixin_41566471/article/details/106319467

最大熵原理

  • 在满足约束条件的模型集合中选取熵最大的模型
  • 离散:

max ⁡ − ∑ x p ( x ) log ⁡ p ( x ) s . t . ∑ x p ( x ) = 1 \begin{aligned} &\max \quad -\sum_x p(x)\log p(x) \\ &s.t. \quad \sum_x p(x)=1 \\ \end{aligned} maxxp(x)logp(x)s.t.xp(x)=1
求解结果是服从均匀分布

  • 连续:连续随机变量的均值为 μ \mu μ,方差为 σ 2 \sigma ^2 σ2,求熵最大对应的概率分布:

max ⁡ L = − ∫ − ∞ ∞ p ( x ) log ⁡ p ( x ) d x s . t . { ∫ − ∞ ∞ p ( x ) d x = 1 ∫ − ∞ ∞ x p ( x ) d x = μ ∫ − ∞ ∞ ( x − μ ) 2 p ( x ) d x = σ 2 \begin{aligned} &\max L=-\int_{-\infty}^{\infty} p(x) \log p(x) d x \\ &s.t. \begin{cases} \int_{-\infty}^{\infty} p(x) d x=1\\ \int_{-\infty}^{\infty} x p(x) d x=\mu\\ \int_{-\infty}^{\infty}(x-\mu)^2 p(x) d x=\sigma^2 \end{cases} \end{aligned} maxL=p(x)logp(x)dxs.t. p(x)dx=1xp(x)dx=μ(xμ)2p(x)dx=σ2
求得模型服从正态分布

最大熵模型的定义

  • 将最大熵原理应用到分类就是最大熵模型。给定训练数据集 T = { ( x 1 , y 1 ) , . . . , ( x N , y N ) } T=\{(x_1,y_1),...,(x_N,y_N)\} T={(x1,y1),...,(xN,yN)},目标就是用最大熵原理选择最好的分类模型。
  • 数学定义为:

image.png

  1. 联合分布 P ( X , Y ) P(X,Y) P(X,Y)和边缘分布 P ( X ) P(X) P(X)的经验分布为:

P ~ ( X = x , Y = y ) = v ( X = x , Y = y ) N \tilde P(X=x,Y=y)=\frac{{\rm v}(X=x,Y=y)}{N} P~(X=x,Y=y)=Nv(X=x,Y=y)
P ~ ( X = x ) = v ( X = x ) N \tilde P(X=x)=\frac{{\rm v}(X=x)}{N} P~(X=x)=Nv(X=x)
v ( X = x , Y = y ) {\rm v}(X=x,Y=y) v(X=x,Y=y)表示训练数据中样本 ( x , y ) (x,y) (x,y)出现的频数。

  1. 特征函数 f ( x , y ) f(x,y) f(x,y)

image.png
注意:特征函数可以是任意实值函数,不一定是二值函数

  1. 特征函数关于经验分布的期望:

E P ~ ( f i ) = ∑ x , y P ~ ( x , y ) f i ( x , y ) E_{\tilde{P}}\left(f_i\right)=\sum_{x, y} \tilde{P}(x, y) f_i(x, y) EP~(fi)=x,yP~(x,y)fi(x,y)

  1. 特征函数关于理论分布的期望:

E P ( f i ) = ∑ x , y P ( x , y ) f i ( x , y ) = ∑ x , y P ( x ) P ( y ∣ x ) f i ( x , y ) = ∑ x , y P ~ ( x ) P ( y ∣ x ) f i ( x , y ) \begin{aligned} E_P\left(f_i\right) &=\sum_{x, y} P(x, y) f_i(x, y) \\ &=\sum_{x, y} P(x) P(y \mid x) f_i(x, y)\\ &=\sum_{x, y}\tilde P(x) P(y \mid x) f_i(x, y) \end{aligned} EP(fi)=x,yP(x,y)fi(x,y)=x,yP(x)P(yx)fi(x,y)=x,yP~(x)P(yx)fi(x,y)

  1. 约束条件:经验分布的期望=理论分布的期望,若有n个特征函数 f i ( x , y ) , i = 1 , . . , n f_i(x,y),i=1,..,n fi(x,y),i=1,..,n,就有n个约束条件。当然,还有一个自然约束条件:概率和为1
  2. 目标:最大熵模型实际上是**「判别方法」**,判别方法就是要得出条件概率分布,就是要得出满足 max ⁡ H ( P ) \max H(P) maxH(P) P ( Y ∣ X ) P(Y|X) P(YX)

H ( P ) = ∑ i = 1 n p ( x i ) H ( Y ∣ X = x i ) = − ∑ i = 1 n p ( x i ) ∑ j = 1 k p ( y j ∣ x i ) log ⁡ p ( y j ∣ x i ) = − ∑ i = 1 n ∑ j = 1 k p ~ ( x i ) p ( y j ∣ x i ) log ⁡ p ( y j ∣ x i ) = − ∑ x , y P ~ ( x ) P ( y ∣ x ) log ⁡ P ( y ∣ x ) \begin{aligned} H(P)&=\sum_{i=1}^n p\left(x_i\right) H\left(Y \mid X=x_i\right)\\ &=-\sum_{i=1}^n p\left(x_i\right)\sum_{j=1}^k p\left(y_j \mid x_i\right) \log p\left(y_j \mid x_i\right) \\ &=-\sum_{i=1}^n \sum_{j=1}^k \tilde{p}\left(x_i\right) p\left(y_j \mid x_i\right) \log p\left(y_j \mid x_i\right)\\ &=-\sum_{x, y} \widetilde{P}(x) P(y \mid x) \log P(y \mid x) \end{aligned} H(P)=i=1np(xi)H(YX=xi)=i=1np(xi)j=1kp(yjxi)logp(yjxi)=i=1nj=1kp~(xi)p(yjxi)logp(yjxi)=x,yP (x)P(yx)logP(yx)

最大熵模型的学习

约束最优化问题

  • 最大熵模型的学习可以形式化为最优化问题:

min ⁡ P ∈ C − H ( P ) = ∑ x , y P ~ ( x ) P ( y ∣ x ) log ⁡ P ( y ∣ x ) s . t . { E P ( f i ) = E P ~ ( f i ) , i = 1 , . . , n ∑ y P ( y ∣ x ) = 1 \begin{aligned} & \min_{P\in C} -H(P)=\sum_{x, y} \widetilde{P}(x) P(y \mid x) \log P(y \mid x)\\ &s.t.\begin{cases} E_P(f_i)=E_{\widetilde P}(f_i),i=1,..,n\\ \sum_y P(y|x)=1 \end{cases} \end{aligned} PCminH(P)=x,yP (x)P(yx)logP(yx)s.t.{EP(fi)=EP (fi),i=1,..,nyP(yx)=1

优化问题的原始问题和对偶问题

image.png

  • f ( x ) , c i ( x ) , h j ( x ) f(x),c_i(x),h_j(x) f(x),ci(x),hj(x) 是连续可微函数。考虑如下有约束的最优化问题:

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{array}{ll} \min _{x \in R^n} & f(x) \\ s . t . & c_i(x) \leq 0, \quad i=1,2, \cdots, k \\ & h_j(x)=0, \quad j=1,2, \cdots, l \end{array} minxRns.t.f(x)ci(x)0,i=1,2,,khj(x)=0,j=1,2,,l

这就是该有约束最优化问题的原始形式。这样一个优化问题其实计算起来非常的复杂,因为有 k+l 个约束条件。

  1. 引进 「广义拉格朗日函数」

L ( x , α , β ) = f ( x ) + ∑ i = 1 k α i c i ( x ) + ∑ j = 1 l β j h j ( x ) α i , β j ≥ 0 L(x, \alpha, \beta)=f(x)+\sum_{i=1}^k \alpha_i c_i(x)+\sum_{j=1}^l \beta_j h_j(x)\quad \alpha_i,\beta_j \ge 0 L(x,α,β)=f(x)+i=1kαici(x)+j=1lβjhj(x)αi,βj0

  1. 「考虑关于x 的函数」:

θ P ( x ) = max ⁡ α , β L ( x , α , β ) = { f ( x ) 有约束 ∞ 无约束 \theta_P(x)=\max _{\alpha, \beta} L(x, \alpha, \beta)=\begin{cases}f(x)\quad 有约束\\ \infty \quad 无约束 \end{cases} θP(x)=α,βmaxL(x,α,β)={f(x)有约束无约束
原始问题的最优值:
p ∗ = min ⁡ x θ P ( x ) = min ⁡ x max ⁡ α , β L ( x , α , β ) p^*=\min _x \theta_P(x)=\min _x \max _{\alpha, \beta} L(x, \alpha, \beta) p=xminθP(x)=xminα,βmaxL(x,α,β)

  1. 「考虑关于 α , β \alpha ,\beta α,β** 的函数」😗*

对偶问题的最优值:
d ∗ = max ⁡ α , β θ D ( α , β ) = max ⁡ α , β min ⁡ x L ( x , α , β ) d^*=\max _{\alpha, \beta} \theta_D(\alpha, \beta)=\max _{\alpha, \beta} \min _x L(x, \alpha, \beta) d=α,βmaxθD(α,β)=α,βmaxxminL(x,α,β)

  1. 原始问题和对偶问题的关系
    1. 广义拉格朗日函数的区间范围

min ⁡ x L ( x , α , β ) ≤ L ( x , α , β ) ≤ max ⁡ α , β L ( x , α , β ) \min _x L(x, \alpha, \beta) \leq L(x, \alpha, \beta) \leq \max _{\alpha, \beta} L(x, \alpha, \beta) xminL(x,α,β)L(x,α,β)α,βmaxL(x,α,β)
d ∗ = max ⁡ α , β : α i ≥ 0 min ⁡ x L ( x , α , β ) ≤ L ( x , α , β ) ≤ min ⁡ x max ⁡ α , β : α i ≥ 0 L ( x , α , β ) = p ∗ d^*=\max _{\alpha, \beta: \alpha_i \geq 0} \min _x L(x, \alpha, \beta) \leq L(x, \alpha, \beta) \leq \min _x \max _{\alpha, \beta: \alpha_i \geq 0} L(x, \alpha, \beta)=p^* d=α,β:αi0maxxminL(x,α,β)L(x,α,β)xminα,β:αi0maxL(x,α,β)=p

  1. 令等号相等需满足的条件

image.png
该条件等价于KKT条件:
∇ x 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 \begin{array}{ll} \nabla_x L\left(x^*, \alpha^*, \beta^*\right)=0 & \\ \alpha_i^* c_i\left(x^*\right)=0, & i=1,2, \cdots, k \\ c_i\left(x^*\right) \leq 0, & i=1,2, \cdots, k \\ \alpha_i^* \geq 0, & i=1,2, \cdots, k \\ h_j(x *)=0, & j=1,2, \cdots, l \end{array} xL(x,α,β)=0αici(x)=0,ci(x)0,αi0,hj(x)=0,i=1,2,,ki=1,2,,ki=1,2,,kj=1,2,,l

证明最大熵模型的原始问题是满足上述条件

image.png
image.png
f ( x ) = − H ( P ) f(x)=-H(P) f(x)=H(P)如果熵函数是严格的凹函数,就要满足:
H ( λ 1 P + λ 2 Q ) > λ 1 H ( P ) + λ 2 H ( Q ) H\left(\lambda_1 P+\lambda_2 Q\right)>\lambda_1 H(P)+\lambda_2 H(Q) H(λ1P+λ2Q)>λ1H(P)+λ2H(Q)
将熵的表达式代入上式:
H ( λ 1 P + λ 2 Q ) = − ∑ i = 1 m ( λ 1 p i + λ 2 q i ) log ⁡ ( λ 1 p i + λ 2 q i ) λ 1 H ( P ) + λ 2 H ( Q ) = − λ 1 ∑ p i log ⁡ p i − λ 2 ∑ q i log ⁡ q i \begin{aligned} &H\left(\lambda_1 P+\lambda_2 Q\right)=-\sum_{i=1}^m\left(\lambda_1 p_i+\lambda_2 q_i\right) \log \left(\lambda_1 p_i+\lambda_2 q_i\right) \\ &\lambda_1 H(P)+\lambda_2 H(Q)=-\lambda_1 \sum p_i \log p_i-\lambda_2 \sum q_i \log q_i \end{aligned} H(λ1P+λ2Q)=i=1m(λ1pi+λ2qi)log(λ1pi+λ2qi)λ1H(P)+λ2H(Q)=λ1pilogpiλ2qilogqi
整理可得我们要证的是下式:
∑ λ 1 p i log ⁡ λ 1 p i + λ 2 q i p i + ∑ λ 2 q i log ⁡ λ 1 p i + λ 2 q i q i < 0 λ 1 + λ 2 = 1 \sum \lambda_1 p_i \log \frac{\lambda_1 p_i+\lambda_2 q_i}{p_i}+\sum \lambda_2 q_i \log \frac{\lambda_1 p_i+\lambda_2 q_i}{q_i}<0 \quad \lambda_1+\lambda_2=1 λ1pilogpiλ1pi+λ2qi+λ2qilogqiλ1pi+λ2qi<0λ1+λ2=1
利用不等式 log ⁡ ( t ) < t − 1 , t ≠ 1 \log (t)<t-1,t \not= 1 log(t)<t1t=1很容易证的。

内部极小化

最大熵模型的原始问题可以通过对偶问题来求解,它的拉格朗日函数为:
L ( P , ω ) = ∑ x , y P ~ ( x ) P ( y ∣ x ) log ⁡ P ( y ∣ x ) + ω 0 ( 1 − ∑ y P ( y ∣ x ) ) + ∑ i = 1 n ω i ( ∑ x , y P ~ ( x , y ) f i ( x , y ) − ∑ x , y P ~ ( x ) P ( y ∣ x ) f i ( x , y ) ) \begin{aligned} L(P, \omega) &=\sum_{x, y} \widetilde{P}(x) P(y \mid x) \log P(y \mid x)+\omega_0\left(1-\sum_y P(y \mid x)\right) \\ &+\sum_{i=1}^n \omega_i\left(\sum_{x, y} \widetilde{P}(x, y) f_i(x, y)-\sum_{x, y} \widetilde{P}(x) P(y \mid x) f_i(x, y)\right) \end{aligned} L(P,ω)=x,yP (x)P(yx)logP(yx)+ω0(1yP(yx))+i=1nωi(x,yP (x,y)fi(x,y)x,yP (x)P(yx)fi(x,y))
其对偶问题为:
d ∗ = max ⁡ ω ≥ 0 θ D ( P ) = max ⁡ ω ≥ 0 min ⁡ P ∈ C L ( P , ω ) d^*=\max _{\omega \geq 0} \theta_D(P)=\max _{\omega \geq 0} \min _{P\in C} L(P,\omega) d=ω0maxθD(P)=ω0maxPCminL(P,ω)
我们先处理 min ⁡ x L ( x , α , β ) \min _x L(x, \alpha, \beta) minxL(x,α,β)这一部分.

  1. ∂ L ( P , ω ) ∂ P ( y ∣ x ) = ∑ x , y P ~ ( x ) [ log ⁡ P ( y ∣ x ) + 1 ] − ω 0 − ∑ i = 1 n ω i ∑ x , y P ~ ( x ) f i ( x , y ) = 0 \frac{\partial L(P, \omega)}{\partial P(y|x)}=\sum_{x, y} \widetilde{P}(x)[\log P(y \mid x)+1]-\omega_0-\sum_{i=1}^n \omega_i \sum_{x, y} \widetilde{P}(x) f_i(x, y)=0 P(yx)L(P,ω)=x,yP (x)[logP(yx)+1]ω0i=1nωix,yP (x)fi(x,y)=0

得到
P ( y ∣ x ) = exp ⁡ ( ∑ i = 1 n w i f i ( x , y ) + ω 0 − 1 ) = exp ⁡ ( ∑ i = 1 n w i f i ( x , y ) ) exp ⁡ ( 1 − ω 0 ) P(y \mid x)=\exp \left(\sum_{i=1}^n w_i f_i(x, y)+\omega_0-1\right)=\frac{\exp \left(\sum_{i=1}^n w_i f_i(x, y)\right)}{\exp \left(1-\omega_0\right)} P(yx)=exp(i=1nwifi(x,y)+ω01)=exp(1ω0)exp(i=1nwifi(x,y))

  1. ∑ y P ( y ∣ x ) = 1 \sum_y P(y|x)=1 yP(yx)=1,得
    exp ⁡ ( 1 − ω 0 ) = ∑ y exp ⁡ ( ∑ i = 1 n w i f i ( x , y ) ) \exp \left(1-\omega_0\right)=\sum_y {\exp \left(\sum_{i=1}^n w_i f_i(x, y)\right)} exp(1ω0)=yexp(i=1nwifi(x,y))
  2. 将2待人1,我们可以得到条件概率分布

P ω ( y ∣ x ) = 1 Z ω ( x ) exp ⁡ ( ∑ i = 1 n w i f i ( x , y ) ) P_{\omega}(y|x)=\frac{1}{Z_{\omega}(x)} \exp \left(\sum_{i=1}^n w_i f_i(x, y)\right) Pω(yx)=Zω(x)1exp(i=1nwifi(x,y))
其中规范化因子
Z ω ( x ) = ∑ y exp ⁡ ( ∑ i = 1 n w i f i ( x , y ) ) Z_\omega(x)=\sum_y \exp \left(\sum_{i=1}^n w_i f_i(x, y)\right) Zω(x)=yexp(i=1nwifi(x,y))

外部极大化

把3中的 P w ( y ∣ x ) P_w(y|x) Pw(yx)带入到拉格朗日函数中,得外部极大化问题为max 对偶函数:
max ⁡ Ψ ( ω ) = ∑ x , y P ~ ( x ) P ω ( y ∣ x ) log ⁡ P ω ( y ∣ x ) + ω 0 ( 1 − ∑ y P ω ( y ∣ x ) ) + ∑ i = 1 n ω i ( ∑ x , y P ~ ( x , y ) f i ( x , y ) − ∑ x , y P ~ ( x ) P ω ( y ∣ x ) f i ( x , y ) ) \begin{aligned} \max \Psi(\omega)=& \sum_{x, y} \widetilde{P}(x) P_\omega(y \mid x) \log P_\omega(y \mid x)+\omega_0\left(1-\sum_y P_\omega(y \mid x)\right) \\ &+\sum_{i=1}^n \omega_i\left(\sum_{x, y} \widetilde{P}(x, y) f_i(x, y)-\sum_{x, y} \widetilde{P}(x) P_\omega(y \mid x) f_i(x, y)\right) \end{aligned} maxΨ(ω)=x,yP (x)Pω(yx)logPω(yx)+ω0(1yPω(yx))+i=1nωi(x,yP (x,y)fi(x,y)x,yP (x)Pω(yx)fi(x,y))

证明对偶函数得极大化=最大熵模型得极大似然估计

将对偶函数继续化简:
∑ y P ω ( y ∣ x ) = 1 \sum_y P_\omega(y \mid x)=1 yPω(yx)=1可知,对偶函数得第二项为0。
Ψ ( ω ) = ∑ x , y P ~ ( x ) P ω ( y ∣ x ) log ⁡ P ω ( y ∣ x ) + ∑ i = 1 n ω i ( ∑ x , y P ~ ( x , y ) f i ( x , y ) − ∑ x , y P ~ ( x ) P ω ( y ∣ x ) f i ( x , y ) ) = ∑ x , y P ~ ( x , y ) ∑ i = 1 n ω i f i ( x , y ) + ∑ x , y P ~ ( x ) P ω ( y ∣ x ) ( log ⁡ P ω ( y ∣ x ) − ∑ i = 1 n ω i f i ( x , y ) ) = ∑ x , y P ~ ( x , y ) ∑ i = 1 n ω i f i ( x , y ) − ∑ x , y P ~ ( x ) P ω ( y ∣ x ) log ⁡ Z ω ( x ) = ∑ x , y P ~ ( x , y ) ∑ i = 1 n ω i f i ( x , y ) − ∑ x P ~ ( x ) log ⁡ Z ω ( x ) \begin{aligned} \Psi(\omega)=& \sum_{x, y} \widetilde{P}(x) P_\omega(y \mid x) \log P_\omega(y \mid x) +\sum_{i=1}^n \omega_i\left(\sum_{x, y} \widetilde{P}(x, y) f_i(x, y)-\sum_{x, y} \widetilde{P}(x) P_\omega(y \mid x) f_i(x, y)\right)\\ =&\sum_{x, y} \widetilde{P}(x, y)\sum_{i=1}^n \omega_i f_i(x, y)+\sum_{x, y} \widetilde{P}(x) P_\omega(y \mid x)\left( \log P_\omega(y \mid x)-\sum_{i=1}^n \omega_if_i(x,y) \right)\\ =&\sum_{x, y} \widetilde{P}(x, y)\sum_{i=1}^n \omega_i f_i(x, y)-\sum_{x, y} \widetilde{P}(x) P_\omega(y \mid x)\log Z_{\omega}(x)\\ =&\sum_{x, y} \widetilde{P}(x, y)\sum_{i=1}^n \omega_i f_i(x, y)-\sum_{x} \widetilde{P}(x) \log Z_{\omega}(x) \end{aligned} Ψ(ω)====x,yP (x)Pω(yx)logPω(yx)+i=1nωi(x,yP (x,y)fi(x,y)x,yP (x)Pω(yx)fi(x,y))x,yP (x,y)i=1nωifi(x,y)+x,yP (x)Pω(yx)(logPω(yx)i=1nωifi(x,y))x,yP (x,y)i=1nωifi(x,y)x,yP (x)Pω(yx)logZω(x)x,yP (x,y)i=1nωifi(x,y)xP (x)logZω(x)
现在看最大熵模型的极大似然估计:
最大熵模型的对数似然函数为:( N ⋅ P ~ ( x , y ) N \cdot \widetilde{P}(x, y) NP (x,y)是数据 ( x , y ) (x,y) (x,y)出现的次数)
L P ~ ( P ω ) = log ⁡ ∏ x , y P ω ( y ∣ x ) N ⋅ P ~ ( x , y ) = N ⋅ log ⁡ ∏ x , y P ω ( y ∣ x ) P ~ ( x , y ) = N ⋅ ∑ x , y P ~ ( x , y ) log ⁡ P ω ( y ∣ x ) = N ⋅ ∑ x , y P ~ ( x , y ) log ⁡ 1 Z ω ( x ) exp ⁡ ( ∑ i = 1 n ω i f i ( x , y ) ) = N [ ∑ x , y P ~ ( x , y ) ( ∑ i = 1 n ω i f i ( x , y ) ) − ∑ x , y P ~ ( x , y ) log ⁡ Z ω ( x ) ] = N [ ∑ x , y P ~ ( x , y ) ( ∑ i = 1 n ω i f i ( x , y ) ) − ∑ x P ~ ( x ) log ⁡ Z ω ( x ) ] \begin{aligned} L_{\widetilde{P}}\left(P_\omega\right) &=\log \prod_{x, y} P_\omega(y \mid x)^{N \cdot \widetilde{P}(x, y)} \\ &=N \cdot \log \prod_{x, y} P_\omega(y \mid x)^{\widetilde{P}(x, y)} \\ &=N \cdot \sum_{x, y} \widetilde{P}(x, y) \log P_\omega(y \mid x)\\ &=N \cdot \sum_{x, y} \widetilde{P}(x, y) \log \frac{1}{Z_\omega(x)} \exp \left(\sum_{i=1}^n \omega_i f_i(x, y)\right) \\ &=N\left[\sum_{x, y} \widetilde{P}(x, y)\left(\sum_{i=1}^n \omega_i f_i(x, y)\right)-\sum_{x, y} \widetilde{P}(x, y) \log Z_\omega(x)\right]\\ &=N\left[\sum_{x, y} \widetilde{P}(x, y)\left(\sum_{i=1}^n \omega_i f_i(x, y)\right)-\sum_{x} \widetilde{P}(x) \log Z_\omega(x)\right] \end{aligned} LP (Pω)=logx,yPω(yx)NP (x,y)=Nlogx,yPω(yx)P (x,y)=Nx,yP (x,y)logPω(yx)=Nx,yP (x,y)logZω(x)1exp(i=1nωifi(x,y))=N[x,yP (x,y)(i=1nωifi(x,y))x,yP (x,y)logZω(x)]=N[x,yP (x,y)(i=1nωifi(x,y))xP (x)logZω(x)]
N是样本数,是固定的,可以省略。这样会发现
Ψ ( ω ) = L P ~ ( P ω ) \Psi(\omega)=L_{\widetilde{P}}\left(P_\omega\right) Ψ(ω)=LP (Pω)

模型学习的最优化算法

我们把求解最大熵模型转化为求解最大化对数似然函数或者最大化对偶函数,那么关键点就转换到找出 拉格朗日乘子 ω \omega ω 上。
max ⁡ Ψ ( ω ) = ∑ x , y P ~ ( x , y ) ∑ i = 1 n ω i f i ( x , y ) − ∑ x P ~ ( x ) log ⁡ Z ω ( x ) \max \Psi(\omega)=\sum_{x, y} \widetilde{P}(x, y)\sum_{i=1}^n \omega_i f_i(x, y)-\sum_{x} \widetilde{P}(x) \log Z_{\omega}(x) maxΨ(ω)=x,yP (x,y)i=1nωifi(x,y)xP (x)logZω(x)
Z ω ( x ) = ∑ y exp ⁡ ( ∑ i = 1 n w i f i ( x , y ) ) Z_\omega(x)=\sum_y \exp \left(\sum_{i=1}^n w_i f_i(x, y)\right) Zω(x)=yexp(i=1nwifi(x,y))

梯度下降法解最大熵模型

目标函数 Q ( ω ) = − Ψ ( ω ) = ∑ x P ~ ( x ) log ⁡ ∑ y exp ⁡ ( ∑ i = 1 n ω i f i ( x , y ) ) − ∑ x , y P ~ ( x , y ) ∑ i = 1 n ω i f i ( x , y ) Q(\omega)=-\Psi(\omega)=\sum_x \widetilde{P}(x) \log \sum_y \exp \left(\sum_{i=1}^n \omega_i f_i(x, y)\right)-\sum_{x, y} \widetilde{P}(x, y) \sum_{i=1}^n \omega_i f_i(x, y) Q(ω)=Ψ(ω)=xP (x)logyexp(i=1nωifi(x,y))x,yP (x,y)i=1nωifi(x,y)

∇ Q ( ω ) = ( ∂ Q ( ω ) ∂ ω 1 , ∂ Q ( ω ) ∂ ω 2 , ⋯   , ∂ Q ( ω ) ∂ ω n ) T \nabla Q(\omega)=\left(\frac{\partial Q(\omega)}{\partial \omega_1}, \frac{\partial Q(\omega)}{\partial \omega_2}, \cdots, \frac{\partial Q(\omega)}{\partial \omega_n}\right)^T Q(ω)=(ω1Q(ω),ω2Q(ω),,ωnQ(ω))T

∂ Q ( ω ) ∂ ω i = ∑ x P ~ ( x ) ∑ y exp ⁡ ( ∑ i = 1 n ω i f i ( x , y ) ) ⋅ f i ( x , y ) ∑ y exp ⁡ ( ∑ i = 1 n ω i f i ( x , y ) ) − ∑ x , y P ~ ( x , y ) f i ( x , y ) = ∑ x P ~ ( x ) P ω ( y ∣ x ) f i ( x , y ) − E P ~ ( f i ) , i = 1 , 2 , ⋯   , n \begin{aligned} \frac{\partial Q(\omega)}{\partial \omega_i}&=\sum_x \widetilde{P}(x) \frac{\sum_y \exp \left(\sum_{i=1}^n \omega_i f_i(x, y)\right) \cdot f_i(x, y)}{\sum_y \exp \left(\sum_{i=1}^n \omega_i f_i(x, y)\right)}-\sum_{x, y} \widetilde{P}(x, y) f_i(x, y)\\ &=\sum_x \widetilde{P}(x) P_\omega(y \mid x) f_i(x, y)-E_{\widetilde{P}}\left(f_i\right), i=1,2, \cdots, n \end{aligned} ωiQ(ω)=xP (x)yexp(i=1nωifi(x,y))yexp(i=1nωifi(x,y))fi(x,y)x,yP (x,y)fi(x,y)=xP (x)Pω(yx)fi(x,y)EP (fi),i=1,2,,n
image.png

拟牛顿法解最大熵模型

牛顿法

先浅谈一下牛顿法:牛顿法最初是为了求解方程的根。

  • 一元

( x ( k ) , g ( x ( k ) ) ) (x^{(k)},g(x^{(k)})) (x(k),g(x(k)))处的切线
y = g ( x ( k ) ) + g ′ ( x ( k ) ) ( x − x ( k ) ) y=g(x^{(k)})+g^{'}(x^{(k)})(x-x^{(k)}) y=g(x(k))+g(x(k))(xx(k))
与 x 轴的交点
y = 0 ⇒ x = x ( k ) − g ( x ( k ) ) g ′ ( x ( k ) ) y=0 \Rightarrow x=x^{(k)}-\frac{g(x^{(k)})}{g^{'}(x^{(k)})} y=0x=x(k)g(x(k))g(x(k))
f ( x ) f(x) f(x)的极值点,就是解 g ( x ) = f ′ ( x ) = 0 g(x)=f^{'}(x)=0 g(x)=f(x)=0的根
x = x ( k ) − g ( x ( k ) ) g ′ ( x ( k ) ) = x ( k ) − f ′ ( x ( k ) ) f ′ ′ ( x ( k ) ) x=x^{(k)}-\frac{g(x^{(k)})}{g^{'}(x^{(k)})}=x^{(k)}-\frac{f^{'}(x^{(k)})}{f^{''}(x^{(k)})} x=x(k)g(x(k))g(x(k))=x(k)f′′(x(k))f(x(k))

  • 多元

x ( k + 1 ) = x ( k ) − H f − 1 ( x ( k ) ) ∇ f ( x ( k ) ) x^{(k+1)}=x^{(k)}-H_f^{-1}\left(x^{(k)}\right) \nabla f\left(x^{(k)}\right) x(k+1)=x(k)Hf1(x(k))f(x(k))
这里的难点就在于要求出 「海森矩阵的逆」,但是如果数值计算量过大,由于是二阶偏导,那在计算过程中就会复杂很多,同时,如果海森矩阵的行列式为零 ,则计算不出逆矩阵,因此这时就需要 「拟牛顿法」出场了~搞清楚了牛顿法的难点,接着我们就看看有什么办法能「替代这个海森矩阵的逆」

拟牛顿条件

牛顿法的实质其实是对它的目标函数 f ( x ) f(x) f(x)进行二阶泰勒展开,这里选择的点就是 x ( k ) x^{(k)} x(k) ,这里另提一句,梯度下降法是通过一阶泰勒展开得到的哦。f对x求导是g
f ( x ) ≈ f ( x ( k ) ) + g k T ( x − x ( k ) ) + 1 2 ( x − x ( k ) ) T H ( x ( k ) ) ( x − x ( k ) ) f(x) \approx f\left(x^{(k)}\right)+g_k^T\left(x-x^{(k)}\right)+\frac{1}{2}\left(x-x^{(k)}\right)^T H\left(x^{(k)}\right)\left(x-x^{(k)}\right) f(x)f(x(k))+gkT(xx(k))+21(xx(k))TH(x(k))(xx(k))
∂ f ( x ) ∂ x ≈ 0 + g ( x ( k ) ) + 1 2 × 2 H ( x ( k ) ) ( x − x ( k ) ) = g k + H k ( x − x ( k ) ) \begin{aligned} \frac{\partial f(x)}{\partial x} & \approx 0+g\left(x^{(k)}\right)+\frac{1}{2} \times 2 H\left(x^{(k)}\right)\left(x-x^{(k)}\right) \\ &=g_k+H_k\left(x-x^{(k)}\right) \end{aligned} xf(x)0+g(x(k))+21×2H(x(k))(xx(k))=gk+Hk(xx(k))
g ( x ( k + 1 ) ) = g k + H k ( x ( k + 1 ) − x ( k ) ) g\left(x^{(k+1)}\right)=g_k+H_k\left(x^{(k+1)}-x^{(k)}\right) g(x(k+1))=gk+Hk(x(k+1)x(k))
「左式」是两个导函数的差值 g k + 1 − g k = y k g_{k+1}-g_{k}=y_k gk+1gk=yk「右式」是两个相邻迭代点之间的差值 x ( k + 1 ) − x ( k ) = δ k x^{(k+1)}-x^{(k)}=\delta_k x(k+1)x(k)=δk
于是可以写成
y k = H k δ k y_k=H_k \delta_k yk=Hkδk
或者
H k − 1 y k = δ k H_k^{-1}y_k=\delta_k Hk1yk=δk
这两个式子成为
拟牛顿条件。

拟牛顿法的合理性

合理性就不证明了,详情看(https://mp.weixin.qq.com/s/NmynLoYV23RCZPXFsYk6PA)

我们知道,拟牛顿法的核心就是为了找到替代品来取代求逆的海森矩阵,那么是不是意味着把这个替代品找到之后,就可以代入牛顿法去计算?
「其实不然!」

拟牛顿法可以看作牛顿法和梯度下降法的结合
{ 牛顿法 : x ( k + 1 ) = x ( k ) − H k − 1 g k 梯度下降法 : x ( k + 1 ) = x ( k ) + λ p k , λ 是步长, p k 是搜索方向 \begin{cases} 牛顿法:x^{(k+1)}=x^{(k)}-H^{-1}_kg_k\\ 梯度下降法:x^{(k+1)}=x^{(k)}+\lambda p_k ,\lambda是步长,p_k是搜索方向 \end{cases} {牛顿法:x(k+1)=x(k)Hk1gk梯度下降法:x(k+1)=x(k)+λpk,λ是步长,pk是搜索方向
拟牛顿法的迭代:步长这样确定 λ = arg min ⁡ λ f ( x ( k ) + λ p k ) \lambda=\argmin_{\lambda}f(x^{(k)}+\lambda p_k) λ=argminλf(x(k)+λpk),搜索方向 p k = − H k − 1 g k p_k=-H^{-1}_kg_k pk=Hk1gk

DFP算法

  • 使用第二个拟牛顿条件就是DFP算法,推导过程如下

image.png
P k y k = δ k ⇒ P k y k = δ k ( δ k T y k ) δ k T y k = δ k δ k T y k δ k T y k ⇒ P k = δ k δ k T δ k T y k P_ky_k=\delta_k \Rightarrow P_ky_k=\frac{\delta_k(\delta_k^Ty_k)}{\delta_k^Ty_k}=\frac{\delta_k\delta_k^Ty_k}{\delta_k^Ty_k} \Rightarrow P_k=\frac{\delta_k\delta_k^T}{\delta_k^Ty_k} Pkyk=δkPkyk=δkTykδk(δkTyk)=δkTykδkδkTykPk=δkTykδkδkT

Q k y k = − G k y k ⇒ Q k y k = − G k y k ( y k T G k y k ) y k T G k y k ⇒ Q k y k = − G k y k y k T G k y k y k T G k y k ⇒ Q k = − G k y k y k T G k y k T G k y k Q_ky_k=-G_ky_k \Rightarrow Q_ky_k=-\frac{G_ky_k(y_k^TG_ky_k)}{y_k^TG_ky_k} \Rightarrow Q_ky_k=-\frac{G_ky_ky_k^TG_ky_k}{y_k^TG_ky_k} \Rightarrow Q_k=-\frac{G_ky_ky_k^TG_k}{y_k^TG_ky_k} Qkyk=GkykQkyk=ykTGkykGkyk(ykTGkyk)Qkyk=ykTGkykGkykykTGkykQk=ykTGkykGkykykTGk
image.png
这里的 G k + 1 G_{k+1} Gk+1就去替代 H k + 1 − 1 H_{k+1}^{-1} Hk+11

DFP求解最大熵模型的完整过程:
image.png
image.png

BFGS算法

  • 使用第一个拟牛顿条件就是BFGS算法,推导过程如下

image.png
这里的 B k + 1 B_{k+1} Bk+1就去替代 H k + 1 H_{k+1} Hk+1

BFGS求解最大熵模型的完整过程:
image.png

Broyden算法

image.png

改进的迭代尺度法IIS

思想:找 δ \delta δ使得 L ( ω + δ ) > L ( ω ) L(\omega+\delta)>L(\omega) L(ω+δ)>L(ω)

在“证明对偶函数得极大化=最大熵模型得极大似然估计”里有:
L ( ω + δ ) − L ( ω ) = ∑ x , y P ~ ( x , y ) log ⁡ P ω + δ ( y ∣ x ) − ∑ x , y P ~ ( x , y ) log ⁡ P ω ( y ∣ x ) = ∑ x , y P ~ ( x , y ) ∑ i = 1 n δ i f i ( x , y ) − ∑ x P ~ ( x ) log ⁡ Z ω + δ ( x ) Z ω ( x ) ≥ ∑ x , y P ~ ( x , y ) ∑ i = 1 n δ i f i ( x , y ) − ∑ x P ~ ( x ) [ 1 − Z ω + δ ( x ) Z ω ( x ) ] = ∑ x , y P ~ ( x , y ) ∑ i = 1 n δ i f i ( x , y ) + 1 − ∑ x P ~ ( x ) Z ω + δ ( x ) Z ω ( x ) \begin{aligned} L(\omega+\delta)-L(\omega)&=\sum_{x, y} \widetilde{P}(x, y) \log P_{\omega+\delta}(y \mid x)-\sum_{x, y} \widetilde{P}(x, y) \log P_\omega(y \mid x)\\ &=\sum_{x, y} \widetilde{P}(x, y) \sum_{i=1}^n \delta_i f_i(x, y)-\sum_x \widetilde{P}(x) \log \frac{Z_{\omega+\delta}(x)}{Z_\omega(x)}\\ &\ge \sum_{x, y} \widetilde{P}(x, y) \sum_{i=1}^n \delta_i f_i(x, y)-\sum_x\widetilde P(x)[1-{\frac{Z_{\omega+\delta}(x)}{Z_\omega(x)}}]\\ &=\sum_{x, y} \widetilde{P}(x, y) \sum_{i=1}^n \delta_i f_i(x, y)+1-\sum_x\widetilde P(x){\frac{Z_{\omega+\delta}(x)}{Z_\omega(x)}} \end{aligned} L(ω+δ)L(ω)=x,yP (x,y)logPω+δ(yx)x,yP (x,y)logPω(yx)=x,yP (x,y)i=1nδifi(x,y)xP (x)logZω(x)Zω+δ(x)x,yP (x,y)i=1nδifi(x,y)xP (x)[1Zω(x)Zω+δ(x)]=x,yP (x,y)i=1nδifi(x,y)+1xP (x)Zω(x)Zω+δ(x)
Z ω + δ ( x ) Z ω ( x ) = ∑ y exp ⁡ ( ∑ i = 1 n ( ω i + δ i ) f i ( x , y ) ) Z ω ( x ) = ∑ y exp ⁡ ( ∑ i = 1 n ω i f i ( x , y ) ) ⋅ exp ⁡ ( ∑ i = 1 n δ i f i ( x , y ) ) Z ω ( x ) = ∑ y P ω ( y ∣ x ) exp ⁡ ∑ i = 1 n δ i f i ( x , y ) \begin{aligned} \frac{Z_{\omega+\delta}(x)}{Z_\omega(x)}&=\frac{\sum_y \exp \left(\sum_{i=1}^n\left(\omega_i+\delta_i\right) f_i(x, y)\right)}{Z_\omega(x)}\\ &=\frac{\sum_y \exp \left(\sum_{i=1}^n \omega_i f_i(x, y)\right) \cdot \exp \left(\sum_{i=1}^n \delta_i f_i(x, y)\right)}{Z_\omega(x)}\\ &=\sum_y P_\omega(y \mid x) \exp \sum_{i=1}^n \delta_i f_i(x, y) \end{aligned} Zω(x)Zω+δ(x)=Zω(x)yexp(i=1n(ωi+δi)fi(x,y))=Zω(x)yexp(i=1nωifi(x,y))exp(i=1nδifi(x,y))=yPω(yx)expi=1nδifi(x,y)

令不等式右边的式子
A ( δ ∣ ω ) = ∑ x , y P ~ ( x , y ) ∑ i = 1 n δ i f i ( x , y ) + 1 − ∑ x P ~ ( x ) ∑ y P ω ( y ∣ x ) exp ⁡ ∑ i = 1 n δ i f i ( x , y ) A(\delta \mid \omega)=\sum_{x, y} \widetilde{P}(x, y) \sum_{i=1}^n \delta_i f_i(x, y)+1-\sum_x \widetilde{P}(x) \sum_y P_\omega(y \mid x) \exp \sum_{i=1}^n \delta_i f_i(x, y) A(δω)=x,yP (x,y)i=1nδifi(x,y)+1xP (x)yPω(yx)expi=1nδifi(x,y)

IIS的真谛

1.找到 A ( δ ∣ ω ) A(\delta \mid \omega) A(δω)** 函数的下界**
引入:
image.png
可以得到如下
{ f i ( x , y ) f # ( x , y ) ≥ 0 ∑ i n f i ( x , y ) f # ( x , y ) = 1 \begin{cases} \begin{gathered} \frac{f_i(x, y)}{f^{\#}(x, y)} \geq 0 \\ \sum_i^n \frac{f_i(x, y)}{f^{\#}(x, y)}=1 \end{gathered} \end{cases} f#(x,y)fi(x,y)0inf#(x,y)fi(x,y)=1
69a24ef1b56c6370a271ef5542323c1.jpg

  1. 说明IIS是合理的

image.png

IIS求解最大熵模型的完整过程

image.png

sklearn没有实现最大熵

  1. 模型实现
import numpy as np
np.random.seed(10)


class MyMaxEntropy(object):
    def __init__(self, lr=0.0001):
        """
        最大熵模型的实现,为了方便理解,尽可能的将参数都存储为字典形式
        :param lr: 学习率,默认值为0.0001

        其他参数:
        :param w: 模型的参数,字典
        :param N: 样本数量
        :param label: 标签空间
        :param hat_p_x: 边缘分布P(X)的经验分布
        :param hat_p_x_y: 联合分布P(X,Y)的经验分布
        :param E_p: 特征函数f(x,y)关于模型P(X|Y)与经验分布hatP(X)的期望值
        :param E_hat_p: 特征函数f(x,y)关于经验分布hatP(X,Y)的期望值
        :param eps: 一个接近于0的正数极小值,这个值放在log的计算中,防止报错
        """
        self.lr = lr
        self.params = {'w': None}

        self.N = None
        self.label = None

        self.hat_p_x = {}
        self.hat_p_x_y = {}

        self.E_p = {}
        self.E_hat_p = {}

        self.eps = np.finfo(np.float32).eps


    def _init_params(self):
        """
        随机初始化模型参数w
        :return:
        """
        w = {}
        for key in self.hat_p_x_y.keys():
            w[key] = np.random.rand()
        self.params['w'] = w

    def _rebuild_X(self, X):
        """
        为了自变量的差异化处理,重新命名自变量
        :param X: 原始自变量
        :return:
        """
        X_result = []
        for x in X:
            X_result.append([y_s + '_' + x_s for x_s, y_s in zip(x, self.X_columns)])
        return X_result

    def _build_mapping(self, X, Y):
        """
        求取经验分布,参照公式(1)(2)
        :param X: 训练样本的输入值
        :param Y: 训练样本的输出值
        :return:
        """
        for x, y in zip(X, Y):
            for x_s in x:
                if x_s in self.hat_p_x.keys():
                    self.hat_p_x[x_s] += 1
                else:
                    self.hat_p_x[x_s] = 1
                if (x_s, y) in self.hat_p_x_y.keys():
                    self.hat_p_x_y[(x_s, y)] += 1
                else:
                    self.hat_p_x_y[(x_s, y)] = 1

        self.hat_p_x = {key: count / self.N for key, count in self.hat_p_x.items()}
        self.hat_p_x_y = {key: count / self.N for key, count in self.hat_p_x_y.items()}

    def _cal_E_hat_p(self):
        """
        计算特征函数f(x,y)关于经验分布hatP(X,Y)的期望值,参照公式(3)
        :return:
        """
        self.E_hat_p = self.hat_p_x_y


    def _cal_E_p(self, X):
        """
        计算特征函数f(x,y)关于模型P(X|Y)与经验分布hatP(X)的期望值,参照公式(4)
        :param X:
        :return:
        """
        for key in self.params['w'].keys():
            self.E_p[key] = 0
        for x in X:
            p_y_x = self._cal_prob(x)
            for x_s in x:
                for (p_y_x_s, y) in p_y_x:
                    if (x_s, y) not in self.E_p.keys():
                        continue
                    self.E_p[(x_s, y)] += (1/self.N) * p_y_x_s

    def _cal_p_y_x(self, x, y):
        """
        计算模型条件概率值,参照公式(9)的指数部分
        :param x: 单个样本的输入值
        :param y: 单个样本的输出值
        :return:
        """

        sum = 0.0
        for x_s in x:
            sum += self.params['w'].get((x_s, y), 0)
        return np.exp(sum), y


    def _cal_prob(self, x):
        """
        计算模型条件概率值,参照公式(9)
        :param x: 单个样本的输入值
        :return:
        """
        p_y_x = [(self._cal_p_y_x(x, y)) for y in self.label]
        sum_y = np.sum([p_y_x_s for p_y_x_s, y in p_y_x])
        return [(p_y_x_s / sum_y, y) for p_y_x_s, y in p_y_x]


    def fit(self, X, X_columns, Y, label, max_iter=20000):
        """
        模型训练入口
        :param X: 训练样本输入值
        :param X_columns: 训练样本的columns
        :param Y: 训练样本的输出值
        :param label: 训练样本的输出空间
        :param max_iter: 最大训练次数
        :return:
        """
        self.N = len(X)
        self.label = label
        self.X_columns = X_columns

        X = self._rebuild_X(X)

        self._build_mapping(X, Y)

        self._cal_E_hat_p()

        self._init_params()

        for iter in range(max_iter):

            self._cal_E_p(X)

            for key in self.params['w'].keys():
                sigma = self.lr * np.log(self.E_hat_p.get(key, self.eps) / self.E_p.get(key, self.eps))
                self.params['w'][key] += sigma

    def predict(self, X):
        """
        预测结果
        :param X: 样本
        :return:
        """
        X = self._rebuild_X(X)
        result_list = []

        for x in X:
            max_result = 0
            y_result = self.label[0]
            p_y_x = self._cal_prob(x)
            for (p_y_x_s, y) in p_y_x:
                if p_y_x_s > max_result:
                    max_result = p_y_x_s
                    y_result = y
            result_list.append((max_result, y_result))
        return result_list

  1. 模型测试
def run_my_model():
    data_set = [['youth', 'no', 'no', '1', 'refuse'],
               ['youth', 'no', 'no', '2', 'refuse'],
               ['youth', 'yes', 'no', '2', 'agree'],
               ['youth', 'yes', 'yes', '1', 'agree'],
               ['youth', 'no', 'no', '1', 'refuse'],
               ['mid', 'no', 'no', '1', 'refuse'],
               ['mid', 'no', 'no', '2', 'refuse'],
               ['mid', 'yes', 'yes', '2', 'agree'],
               ['mid', 'no', 'yes', '3', 'agree'],
               ['mid', 'no', 'yes', '3', 'agree'],
               ['elder', 'no', 'yes', '3', 'agree'],
               ['elder', 'no', 'yes', '2', 'agree'],
               ['elder', 'yes', 'no', '2', 'agree'],
               ['elder', 'yes', 'no', '3', 'agree'],
               ['elder', 'no', 'no', '1', 'refuse'],
               ]
    columns = ['age', 'working', 'house', 'credit_situation', 'label']
    X = [i[:-1] for i in data_set]
    X_columns = columns[:-1]
    Y = [i[-1] for i in data_set]
    print(X)
    print(Y)

    my = MyMaxEntropy()
    train_X = X[:12]
    test_X = X[12:]
    train_Y = Y[:12]
    test_Y = Y[12:]
    my.fit(train_X, X_columns, train_Y, label=['refuse', 'agree'])

    print(my.params)

    pred_Y= my.predict(test_X)
    print('result: ')
    print('test: ', test_Y)
    print('pred: ', pred_Y)

  • 0
    点赞
  • 2
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

weixin_961876584

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值