超详尽的变分推断算法教程及例子

k1、变分推断的算法

在推导之前先陈述一下各个变量的含义

X:观测变量,Z:隐变量+参数

根据我们的贝叶斯公式,我们所要求的后验分布为:
p ( Z ∣ X ) = p ( X , Z ) p ( X ) p(Z|X) = \frac {p(X,Z)}{p(X)} p(ZX)=p(X)p(X,Z)

p ( Z ∣ X ) = p ( X , Z ) ∫ Z p ( X ∣ Z ) p ( Z ) d Z p(Z|X) = \frac {p(X,Z)}{\int_Z p(X|Z)p(Z)dZ} p(ZX)=Zp(XZ)p(Z)dZp(X,Z)

很多时候 p ( Z ∣ X ) p(Z|X) p(ZX)是很不好求的原因是当Z的维度很大时分母进行积分就会变的很复杂求不出解析解因此我们就需要通过用一个 q ( Z ) q(Z) q(Z)的分布来近似真实的 p ( Z ) p(Z) p(Z)分布,我们可以将上述的式子进行化简,引入 q ( z ) q(z) q(z)

进行简单的变换我们可以得到
p ( X ) = p ( X , Z ) p ( Z ∣ X ) p(X) = \frac {p(X,Z)}{p(Z|X)} p(X)=p(ZX)p(X,Z)
两边同时取对数可以得到:
l o g P ( X ) = l o g p ( X , Z ) p ( Z ∣ X ) logP(X) = log\frac {p(X,Z)}{p(Z|X)} logP(X)=logp(ZX)p(X,Z)

= l o g p ( X , Z ) − l o g p ( Z ∣ X ) =log{p(X,Z)}-log{p(Z|X)} =logp(X,Z)logp(ZX)

= l o g p ( X , Z ) q ( Z ) − l o g p ( Z ∣ X ) q ( Z ) =log\frac {p(X,Z)}{q(Z)}-log\frac {p(Z|X)}{q(Z)} =logq(Z)p(X,Z)logq(Z)p(ZX)

= l o g p ( X , Z ) − l o g q ( Z ) − l o g p ( Z ∣ X ) q ( Z ) =log{p(X,Z)}-log{q(Z)}-log\frac {p(Z|X)}{q(Z)} =logp(X,Z)logq(Z)logq(Z)p(ZX)

最后一行的 q ( X ) q(X) q(X)是我们需要用来拟合 p ( X ) p(X) p(X)的一个分布

1.1公式化简

我们对上述公式同时左右两边同时求q(Z)期望可以得到

左边:

对于左边的式子 l o g p ( X ) logp(X) logp(X)求q(Z)的期望得:
E ( l o g p ( X ) ) = ∫ Z l o g p ( X ) q ( Z ) d Z = l o g p ( X ) ∫ Z q ( Z ) d Z E(logp(X))=\int_Z logp(X)q(Z)dZ=logp(X)\int_Z q(Z)dZ E(logp(X))=Zlogp(X)q(Z)dZ=logp(X)Zq(Z)dZ
由于 ∫ Z q ( Z ) d Z = 1 \int_Z q(Z)dZ=1 Zq(Z)dZ=1故最终的 E ( l o g p ( x ) ) = l o g p ( X ) E(logp(x))=logp(X) E(logp(x))=logp(X)

右边:

l o g p ( X , Z ) − l o g q ( Z ) − l o g p ( Z ∣ X ) q ( Z ) = ∫ Z l o g p ( X , Z ) q ( Z ) d Z − ∫ Z l o g q ( Z ) q ( Z ) d Z ⏟ E L B O + ∫ Z − l o g p ( Z ∣ X ) q ( Z ) q ( Z ) d Z ⏟ K L log{p(X,Z)}-log{q(Z)}-log\frac {p(Z|X)}{q(Z)}=\underbrace {\int_Z log{p(X,Z)}q(Z)dZ-\int_Z logq(Z)q(Z)dZ}_{ELBO}+\underbrace{\int_Z -log\frac {p(Z|X)}{q(Z)}q(Z)dZ}_{KL} logp(X,Z)logq(Z)logq(Z)p(ZX)=ELBO Zlogp(X,Z)q(Z)dZZlogq(Z)q(Z)dZ+KL Zlogq(Z)p(ZX)q(Z)dZ

= ∫ Z l o g p ( X , Z ) q ( Z ) q ( Z ) d Z ⏟ E L B O + ∫ Z − l o g p ( Z ∣ X ) q ( Z ) q ( Z ) d Z ⏟ K L =\underbrace{\int_Z log\frac {p(X,Z)}{q(Z)}q(Z)dZ}_{ELBO}+\underbrace{\int_Z -log\frac {p(Z|X)}{q(Z)}q(Z)dZ}_{KL} =ELBO Zlogq(Z)p(X,Z)q(Z)dZ+KL Zlogq(Z)p(ZX)q(Z)dZ

最终化简的形式如下:
l o g p ( X ) = ∫ Z l o g p ( X , Z ) q ( Z ) q ( Z ) d Z ⏟ E L B O + ∫ Z − l o g p ( Z ∣ X ) q ( Z ) q ( Z ) d Z ⏟ K L logp(X)=\underbrace{\int_Z log\frac {p(X,Z)}{q(Z)}q(Z)dZ}_{ELBO}+\underbrace{\int_Z -log\frac {p(Z|X)}{q(Z)}q(Z)dZ}_{KL} logp(X=ELBO Zlogq(Z)p(X,Z)q(Z)dZ+KL Zlogq(Z)p(ZX)q(Z)dZ

其中 E L B O ELBO ELBO也是被记为 L ( q ) L(q) L(q)

而KL散度可以衡量两个分布 p p p q q q之间的距离

在这里插入图片描述

此公式的含义也可以由上图所表示出来

由于我们求不出 p ( Z ∣ X ) p(Z|X) p(ZX),我们的目的是寻找一个 q ( Z ) q(Z) q(Z)使得 p ( Z ∣ X ) p(Z|X) p(ZX)近似 q ( Z ) q(Z) q(Z)也就是 K L ( q ∣ ∣ p ( z ∣ x ) ) KL(q||p(z|x)) KL(qp(zx))越小越好(p和q越相近越好,其中为什么是近似 p ( z ∣ x ) p(z|x) p(zx)而不是 p ( z ) p(z) p(z)的原因是 p ( z ) p(z) p(z)很多时候并不能求出来具体的值,另一方面就是我们可以利用x的值来求后验概率缩小概率的范围)。并且,$p(X) 是 个 定 值 ( 原 因 是 P ( X ) 代 表 真 实 的 分 布 , 真 实 分 布 是 存 在 的 也 是 一 个 具 体 的 值 , 只 不 过 是 我 们 不 知 道 它 具 体 的 值 是 多 少 , 因 此 对 应 的 X 的 概 率 是 一 个 具 体 的 值 ) , 那 么 我 们 的 目 标 变 成 了 是个定值(原因是P(X)代表真实的分布,真实分布是存在的也是一个具体的值,只不过是我们不知道它具体的值是多少,因此对应的X的概率是一个具体的值),那么我们的目标变成了 P(X)X argmax_{q(z)}L(q)$。那么,我们理一下思路,

我们想要求得一个 q ( Z ) ≈ p ( Z ∣ X ) q(Z)≈p(Z|X) q(Z)p(ZX)。也就是
q ( Z ) = a r g m a x q ( z ) L ( q ) q(Z)=argmax_{q(z)}L(q) q(Z)=argmaxq(z)L(q)

1.2模型求解(经典平均场理论)

那么我们如何来求解这个问题呢?我们使用到统计物理中的一种方法,就是平均场理论 (mean field theory)。也就是假设变分后验分式是一种完全可分解的分布:

q ( z ) = ∏ i = 1 M q i ( z i ) q(z)=\prod_{i=1}^Mq_i(z_i) q(z)=i=1Mqi(zi)
也就是认为每个 z i z_i zi都是相互独立的(但是这种方法的条件性很强,造成的问题就是对于一些比较复杂的问题拟合的 p ( z ) p(z) p(z)会有比较大的bias,好处就是将高维的参数转化成了低纬度,极大简化计算)我们将 q ( z ) q(z) q(z)带入到 L ( q ) L(q) L(q)

通过平均场理论我们就把原先高维的隐变量转化成低维的独立的隐变量相乘,极大简化的积分的复杂度
L ( q ) = ∫ Z l o g p ( X , Z ) q ( Z ) d Z − ∫ Z l o g q ( Z ) q ( Z ) d Z L(q)=\int_Z log{p(X,Z)}q(Z)dZ-\int_Z logq(Z)q(Z)dZ L(q)=Zlogp(X,Z)q(Z)dZZlogq(Z)q(Z)dZ

= ∫ Z l o g p ( X , Z ) ∏ i = 1 M q i ( z i ) d Z − ∫ Z ∑ i = 1 M l o g q ( Z ) ∏ i = 1 M q i ( z i ) d Z =\int_Z log{p(X,Z)}\prod_{i=1}^Mq_i(z_i)dZ-\int_Z \sum _{i=1}^Mlogq(Z)\prod_{i=1}^Mq_i(z_i)dZ =Zlogp(X,Z)i=1Mqi(zi)dZZi=1Mlogq(Z)i=1Mqi(zi)dZ

然后我们将上述式子分成两个部分

part 1:

∫ Z l o g p ( X , Z ) q ( Z ) d Z = ∫ Z l o g p ( X , Z ) ∏ i = 1 M q i ( z i ) d Z \int_Z log{p(X,Z)}q(Z)dZ=\int_Z log{p(X,Z)}\prod_{i=1}^Mq_i(z_i)dZ Zlogp(X,Z)q(Z)dZ=Zlogp(X,Z)i=1Mqi(zi)dZ

我们将整个集合的Z拆解成一个个z进行积分
= ∫ z 1 ∫ z 2 . . . ∫ z M q i ( z i ) l o g ( p ( X , Z ) ) d z 1 d z 2 . . . d z M =\int_{z_1}\int_{z_2}...\int_{z_M}q_i(z_i)log(p(X,Z))dz_1dz_2...dz_M =z1z2...zMqi(zi)log(p(X,Z))dz1dz2...dzM
关键的地方来了,我们这里除了第 Z j Z_j Zj个隐变量讲其他项的都积分掉,这样就可以留下含有 Z j Z_j Zj项的式子,这样我们后面就可以通过CAVI(坐标上升法)来对每个参数进行迭代更新,
= ∫ z j q j ( z j ) [ ∫ z 1 ∫ z 2 . . . ∫ z M q i ( z i ) l o g ( p ( X , Z ) ) d z 1 d z 2 . . . d z M ] d j =\int_{z_j}q_j(z_j)[\int_{z_1}\int_{z_2}...\int_{z_M}q_i(z_i)log(p(X,Z))dz_1dz_2...dz_M]d_j =zjqj(zj)[z1z2...zMqi(zi)log(p(X,Z))dz1dz2...dzM]dj

= ∫ z j q j ( z j ) E ∏ i ≠ j M q i ( z i ) [ l o g p ( X , Z ) ] d z j =\int_{z_j}q_j(z_j)E_{\prod_{i\neq j}^Mq_i(z_i)}[logp(X,Z)]dz_j =zjqj(zj)Ei=jMqi(zi)[logp(X,Z)]dzj

part 2:

∫ Z ∑ i = 1 M l o g q ( Z ) ∏ i = 1 M q i ( z i ) d Z = ∫ Z ∏ i = 1 M q i ( z i ) d Z [ l o g q 1 ( z 1 ) + l o g q 2 ( z 2 ) + . . . + l o g q M ( z M ) ] d Z \int_Z \sum _{i=1}^Mlogq(Z)\prod_{i=1}^Mq_i(z_i)dZ=\int_Z \prod_{i=1}^Mq_i(z_i)dZ[logq_1(z_1)+logq_2(z_2)+...+logq_M(z_M)]dZ Zi=1Mlogq(Z)i=1Mqi(zi)dZ=Zi=1Mqi(zi)dZ[logq1(z1)+logq2(z2)+...+logqM(zM)]dZ

我们如何计算呢?现在我们通过简单的式子来化简出通向公式,现在设定 M = 2 M=2 M=2
∫ z 1 ∫ z 2 [ l o g q 1 ( z 1 ) + l o g q 2 ( z 2 ) ] q 1 ( z 1 ) q 2 ( z 2 ) d z 1 d z 2 = ∫ z 1 ∫ z 2 q 1 ( z 1 ) q 2 ( z 2 ) l o g q 1 ( z 1 ) d z 1 d z 2 + ∫ z 1 ∫ z 2 q 1 ( z 1 ) q 2 ( z 2 ) l o g q 2 ( z 2 ) d z 1 d z 2 \int_{z_1}\int_{z_2}[logq_1(z_1)+logq_2(z_2)]q_1(z_1)q_2(z_2)dz_1dz_2=\int_{z_1}\int_{z_2}q_1(z_1)q_2(z_2)logq_1(z_1)dz_1dz_2+\int_{z_1}\int_{z_2}q_1(z_1)q_2(z_2)logq_2(z_2)dz_1dz_2 z1z2[logq1(z1)+logq2(z2)]q1(z1)q2(z2)dz1dz2=z1z2q1(z1)q2(z2)logq1(z1)dz1dz2+z1z2q1(z1)q2(z2)logq2(z2)dz1dz2

= ∫ z 1 q 1 ( z 1 ) l o g q 1 ( z 1 ) ∫ z 2 q 2 ( z 2 ) d z 2 ⏟ 1 d z 1 + ∫ z 2 q 2 ( z 2 ) l o g q 2 ( z 2 ) ∫ z 1 q 1 ( z 1 ) d z 1 ⏟ 1 d z 2 =\int_{z_1}q_1(z_1)logq_1(z_1)\underbrace{\int_{z_2}q_2(z_2)dz_2}_1dz_1+\int_{z_2}q_2(z_2)logq_2(z_2)\underbrace{\int_{z_1}q_1(z_1)dz_1}_1dz_2 =z1q1(z1)logq1(z1)1 z2q2(z2)dz2dz1+z2q2(z2)logq2(z2)1 z1q1(z1)dz1dz2

= ∑ i = 1 2 ∫ z i q i ( z i ) l o g q i ( z i ) d z i =\sum_{i=1}^2\int_{z_i}q_i(z_i)logq_i(z_i)dz_i =i=12ziqi(zi)logqi(zi)dzi

故可以求得通向公式
∫ Z ∑ i = 1 M l o g q ( Z ) ∏ i = 1 M q i ( z i ) d Z = ∑ i = 1 M ∫ z i q i ( z i ) l o g q i ( z i ) d z i \int_Z \sum _{i=1}^Mlogq(Z)\prod_{i=1}^Mq_i(z_i)dZ=\sum_{i=1}^M\int_{z_i}q_i(z_i)logq_i(z_i)dz_i Zi=1Mlogq(Z)i=1Mqi(zi)dZ=i=1Mziqi(zi)logqi(zi)dzi
因为我们仅仅只关注第 j 项因为第j项是我们需要进行迭代的一项,所以除了j项其余项都为常数
p a r t 2 = ∫ z j q i ( z i ) l o g q i ( z i ) d z j + C part2=\int_{z_j}q_i(z_i)logq_i(z_i)dz_j+C part2=zjqi(zi)logqi(zi)dzj+C

为了近一步计算我们将上述 p a r t 1 part 1 part1式子 E ∏ i ≠ j M q i ( z i ) [ l o g p ( X , Z ) ] E_{\prod_{i\neq j}^Mq_i(z_i)}[logp(X,Z)] Ei=jMqi(zi)[logp(X,Z)]表示为
E ∏ i ≠ j M q i ( z i ) [ l o g p ( X , Z ) ] = l o g p ^ ( X , z j ) E_{\prod_{i\neq j}^Mq_i(z_i)}[logp(X,Z)]=log\hat p(X,z_j) Ei=jMqi(zi)[logp(X,Z)]=logp^(X,zj)
为什么要这样做呢,目的是为了可以构造出 L ( q ) = − K L L(q)=-KL L(q)=KL后面的式子,求出我们每个 Z j Z_j Zj分布需要逼近的真实的分布,我们可以对上述式子求逆得到伪真实分布 p ^ ( X , z j ) \hat p(X,z_j) p^(X,zj)为什么可以叫做伪真实分布因为它们之间只是相差了一个常数项
p ^ ( X , z j ) = e x p E ∏ i ≠ j M q i ( z i ) [ l o g p ( X , Z ) ] d z j \hat p(X,z_j)=exp^{E_{\prod_{i\neq j}^Mq_i(z_i)}[logp(X,Z)]dz_j} p^(X,zj)=expEi=jMqi(zi)[logp(X,Z)]dzj
那么 p a r t 1 part 1 part1的式子就可以写做
= ∫ z j q j ( z j ) l o g p ^ ( X , z j ) d z j =\int_{z_j}q_j(z_j)log\hat p(X,z_j)dz_j =zjqj(zj)logp^(X,zj)dzj
这样我们就可以求得最终的 L ( q ) L(q) L(q)的表达式
L ( q j ) = p a r t 1 − p a r t 2 = ∫ z j q j ( z j ) l o g p ^ ( X , z j ) d z j − ∫ z j q i ( z j ) l o g q i ( z j ) d z j + C L(q_j)=part1-part2=\int_{z_j}q_j(z_j)log\hat p(X,z_j)dz_j-\int_{z_j}q_i(z_j)logq_i(z_j)dz_j+C L(qj)=part1part2=zjqj(zj)logp^(X,zj)dzjzjqi(zj)logqi(zj)dzj+C

= ∫ z j q j ( z j ) l o g p ^ ( X , z j ) q i ( z j ) d z j + C = − K L =\int_{z_j}q_j(z_j)log\frac{\hat p(X,z_j)}{q_i(z_j)}dz_j+C=-KL =zjqj(zj)logqi(zj)p^(X,zj)dzj+C=KL

= ∫ z j q j ( z j ) l o g p ^ ( X , z j ) d z j + C =\int_{z_j}q_j(z_j)log\hat p(X,z_j)dz_j+C =zjqj(zj)logp^(X,zj)dzj+C

因此最后我们要求解 a r g m a x q ( z ) L ( q j ) argmax_{q(z)}L(q_j) argmaxq(z)L(qj)就是求解等价于 a r g m i n q ( z ) K L ( q j ∣ ∣ p ^ ( x , z j ) ) argmin_{q(z)}KL(q_j||\hat p(x,z_j)) argminq(z)KL(qjp^(x,zj)) ,有kl散度的定义我们可以知道只有当 q j q_j qj等于 p ^ ( x , z j ) \hat p(x,z_j) p^(x,zj)时可以达到最优解,对此我们就可以得到每个隐变量的迭代更新的式子
q i ( z i ) = e x p E ∏ i ≠ j M q i ( z i ) [ l o g p ^ ( X , Z ) ] d z j q_i(z_i) = exp^{E_{\prod_{i\neq j}^Mq_i(z_i)}[log\hat p(X,Z)]dz_j} qi(zi)=expEi=jMqi(zi)[logp^(X,Z)]dzj

l o g q i ( z i ) = E ∏ i ≠ j M q i ( z i ) [ l o g p ^ ( X , Z ) ] d z j logq_i(z_i) =E_{\prod_{i\neq j}^Mq_i(z_i)}[log\hat p(X,Z)]dz_j logqi(zi)=Ei=jMqi(zi)[logp^(X,Z)]dzj

下面以一个业务例子来实践这个变分推断的方法。

1.3一个业务的例子

首先,于是我们根据业务知识,知道用户产生时长增量记录的过程为:

假设在我们的app中总用户数为k,并且在app的数据库中记录了用户每日使用时长的增长量,假设我们没有任何关于用户的唯一id,能观测到的只有增量值: [-6.04231708, -1.64784446, …, 1.63898137, -4.29063685, -0.87892713] ,我们需要做的是为每个用户赋予一个用户id,并且在未来时刻给定某个用户使用时长增量,将这个时长增量归属到其用户id上,如此便可以建立每个用户的使用时长记录,以便在商业上分析用户行为。

  • 某个用户 c k c_k ck登录系统
  • 用户产生一条使用时长增量记录 x i x_i xi

然后,我们将这个过程建模为数学问题:假设登录到系统上的用户为 k ′ k' k的真实分布是均匀分布的,并且概率都为 1 k \frac 1k k1并且用户的 k ′ k' k的增量为一个随机变量,其符合均值为 u k u_k uk标准差为1的分布特别注意的是 σ \sigma σ是一个人工设定的超参数,为领域专家根据先验知识调整的。数学化的表述如下:

μ k \mu_k μk代表的是我的观测数据 x i x_i xi所属的分布的均值
μ k ∼ N ( 0 , σ 2 )         k = 1 , . . . . K \mu_k \sim N(0,\sigma^2)\ \ \ \ \ \ \ k=1,....K μkN(0,σ2)       k=1,....K
我们可以写出 u k u_k uk的概率密度函数
p ( u k ) = 1 2 π σ 2 e − u k 2 2 σ 2 p(u_k)=\frac {1}{\sqrt {2\pi\sigma^2}}e^{-\frac {u_k^2}{2\sigma^2}} p(uk)=2πσ2 1e2σ2uk2
c i c_i ci服从类别分布
c i ∼ C a t e g o r i c a l ( 1 k , 1 k , . . . 1 k )          i = 1 , . . . . . , n c_i\sim Categorical(\frac 1k,\frac 1k,...\frac 1k)\ \ \ \ \ \ \ \ i=1,.....,n ciCategorical(k1,k1,...k1)        i=1,.....,n
p ( c i ) p(c_i) p(ci)的概率密度函数如下:
p ( c i ) = 1 k p(c_i) = \frac {1}{k} p(ci)=k1
c i c_i ci代表的是第i条观测数据 x i x_i xi是来自于哪个分布,例如:当K=6时 c i = [ 0 , 0 , 1 , 0 , 0 , 0 ] ci=[0,0,1,0,0,0] ci=[0,0,1,0,0,0]的向量,代表第i条数据 x i x_i xi属于第3个分布
x i ∣ c i ∼ N ( c i , μ , 1 )          i = 1 , . . . . , n x_i|c_i \sim N(c_i,\mu,1)\ \ \ \ \ \ \ \ i=1,....,n xiciN(ci,μ,1)        i=1,....,n

值得注意的是由于我们的假设 p ( μ k ) p(\mu_k) p(μk) p ( c i ) p(c_i) p(ci)是相互独立的先验概率,后面我们需要通过一个我们定义的 q ( z ) q(z) q(z)分布去拟合这些真实的分布 p ( z ∣ x ) p(z|x) p(zx)注意这里是给定x条件的z的分布,这样的目的是我们可以把原先的先验分布转化成后验分布,那么去拟合后验分布有什么好处呢,就举一个简单的例子,对于判断给定的变量x是来自哪个分布的 u k u_k uk如果按照先验分布去拟合,最后得到的结果就是对于每个变量x来自于每个 u k u_k uk的概率都是一样的,这显然是毫无意义的,我们的目的就是为了让模型可以帮助我们把变量x服从哪个分布来区分出来,那我们就需要依靠现有的数据集X来作为我的条件求解z的后验分布,不断更新z的值从而能够使模型能够区分给定x是来自于哪个分布

根据上述我们设定的这些数据的真实分布我们可以利用这些分布来生成一些伪造的真实数据,我们可以用python代码写出整个数据产生的过程:

def gen_data(sigma, k, n):
    #获取u
    u = np.random.normal(0, sigma, k)
    x = []
    c = []
    for i in range(n):
        ci = random.choice(range(k))
        c.append(ci)
        u = u_k[ci]
        x.append(np.random.normal(u, 1))
    return x, u_k, c

x, u, c = gen_data(sigma,k,n)

我们使用这个函数产生 ( σ \sigma σ=10, k=2)实验室数据,采样得到的

x:[3.91063798840766, 17.707494246717758, 5.1750776782449845, 16.97945369198913, 5.386555860818202, 16.627103391787628, 5.043095807216062, 17.35238536362031, 18.639211648953882, 15.400741385911331]

u u u:[16.93029991 4.53146332]
x所属分布class:[1, 0, 1, 0, 1, 0, 1, 0, 0, 0](因为这里k为2故列表的范围只能是0和1,分别代表两个分布,例如class[0]代表第一个数据x[0]服从分布1)

以下的通用公式都会通过这个例子来进行解释

数据x分布如下图:

在这里插入图片描述

于是我们的目标是根据数据x,去计算得出u,然后未来有数据x′到来时查看数据与 u i ∈ k u_{i\in k} uik的距离,选择距离最近的$ i=argmin_{i} dis(u_{i \in k}, x’)$赋予该条数据作为其用户id即可。

接下来我们来求我们需要求的隐变量 u k , c i u_k,c_i uk,ci所服从的分布对应的参数的值,这样当我们求解得出这些分布的参数后,我们就可以通过这些分布参数来求解这些隐变量最有可能的取值。从而推测出这些隐变量的值,对于上述隐变量我们需要求[ φ 1 , 0 , . . . , φ 10 , 0 , φ 1 , 1 , . . . , φ 10 , 1 , m 0 , s 0 , m 1 , s 1 \varphi_{1,0},...,\varphi_{10,0},\varphi_{1,1},...,\varphi_{10,1},m_0,s_0,m_1,s_1 φ1,0,...,φ10,0,φ1,1,...,φ10,1,m0,s0,m1,s1]总共24个参数

为了求解第三节中的问题,根据平均场变分推断的思路,我们先假设隐变量 u k u_k uk的分布满足均值为 m k m_k mk方差为 s k s_k sk而隐变量 c i c_i ci由参数 φ i \varphi_i φi决定, 参数也是k维向量,每一维度指定了 c i c_i ci对应维度为1的概率,即:
q ( c i = k ; φ i ) = φ i k q(c_i=k;\varphi_i)=\varphi_{ik} q(ci=k;φi)=φik

φ i k \varphi_{ik} φik代表的是属于该分布的概率,对应的因为我的c是满足分类分布 c i k c_{ik} cik的取值范围只能是0或1, , q ( c i = 2 ; φ i ) = φ i 2 , q ( c i = 1 ; φ i ) = φ i 1 ,q(c_i=2;\varphi_i)=\varphi_{i2},q(c_i=1;\varphi_i)=\varphi_{i1} q(ci=2;φi)=φi2q(ci=1;φi)=φi1

先求解 q ( c i ; φ i ) q(c_i;\varphi_i) q(ci;φi)

有上述证明的 l o g q j ∗ ( z j ) ∝ E ∏ i ≠ j M q i ( z i ) [ l o g p ( z , x ) ] logq^*_j(z_j) \propto E_{\prod_{i\neq j}^Mq_i(z_i)}[logp(z,x)] logqj(zj)Ei=jMqi(zi)[logp(z,x)](只有通过使用平均随机场时才可以用这个公式,

这里我们回到本例例子,求 q ∗ ( c i ; φ i ) q^*(c_i;\varphi_i) q(ci;φi),带入公式得
l o g q ∗ ( c i ; φ i ) ∝ E ∏ j ≠ i M q j ( z j ) [ l o g p ( c , u , x ) ] logq^*(c_i;\varphi_i)\propto E_{\prod_{j\neq i}^Mq_j(z_j)}[log p(c,u,x)] logq(ci;φi)Ej=iMqj(zj)[logp(c,u,x)]
这里我们的 q j ( z j ) q_j(z_j) qj(zj)代表 q 1 ( c 1 ) , q 2 ( c 2 ) , . . . . q 10 ( c 10 ) , q 11 ( u 1 ) q 12 ( u 2 ) q_1(c_1),q_2(c_2),....q_{10}(c_{10}),q_{11}(u_1)q_{12}(u_2) q1(c1),q2(c2),....q10(c10),q11(u1)q12(u2)共计12个z
E ∏ j ≠ i M q j ( z j ) [ l o g p ( c , u , x ) ] = E q − c i [ l o g [ p ( μ ) p ( c ) p ( x ∣ μ , c ) ] ] E_{\prod_{j\neq i}^Mq_j(z_j)}[log p(c,u,x)] = E_{q-c_i}[log[p(\mu)p(c)p(x|\mu,c)]] Ej=iMqj(zj)[logp(c,u,x)]=Eqci[log[p(μ)p(c)p(xμ,c)]]

= E q − c i [ l o g [ p ( μ ) p ( c i ) p ( c − i ) p ( x i ∣ μ , c i ) p ( x − i ∣ μ , c − i ) ] ] =E_{q-c_i}[log[p(\mu)p(c_i)p(c_{-i})p(x_i|\mu,c_i)p(x_{-i}|\mu,c_{-i})]] =Eqci[log[p(μ)p(ci)p(ci)p(xiμ,ci)p(xiμ,ci)]]

= E q − c i [ l o g p ( x i ∣ μ , c i ) ] + E q − c i [ l o g p ( μ ) ] + E q − c i [ l o g p ( c i ) ] + E q − c i [ l o g p ( c − i ) p ( x − i ∣ μ , c − i ) ] =E_{q-c_i}[logp(x_i|\mu,c_i)]+E_{q-c_i}[logp(\mu)]+E_{q-c_i}[logp(c_i)]+E_{q-c_i}[logp(c_{-i})p(x_{-i}|\mu,c_{-i})] =Eqci[logp(xiμ,ci)]+Eqci[logp(μ)]+Eqci[logp(ci)]+Eqci[logp(ci)p(xiμ,ci)]

= E q − c i [ l o g p ( x i ∣ μ , c i ) ] + E q − c i [ l o g p ( μ ) ] + l o g p ( c i ) + E q − c i [ l o g p ( c − i ) p ( x − i ∣ μ , c − i ) ] =E_{q-c_i}[logp(x_i|\mu,c_i)]+E_{q-c_i}[logp(\mu)]+logp(c_i)+E_{q-c_i}[logp(c_{-i})p(x_{-i}|\mu,c_{-i})] =Eqci[logp(xiμ,ci)]+Eqci[logp(μ)]+logp(ci)+Eqci[logp(ci)p(xiμ,ci)]

从上述推导中,仅仅只有 E q − c i [ l o g p ( x i ∣ μ , c i ) ] E_{q-c_i}[logp(x_i|\mu,c_i)] Eqci[logp(xiμ,ci)] l o g p ( c i ) logp(c_i) logp(ci)是依赖于 c i c_i ci的项,其余项都可以提出为常数项,因此我们就可以更新我们的正比公式
l o g q ∗ ( c i ; φ i ) ∝ E q − c i [ l o g p ( x i ∣ μ , c i ) ] + l o g p ( c i ) logq^*(c_i;\varphi_i)\propto E_{q-c_i}[logp(x_i|\mu,c_i)]+logp(c_i) logq(ci;φi)Eqci[logp(xiμ,ci)]+logp(ci)
通过先验,我们可以得到
l o g p ( c i ) = − l o g K logp(c_i) = -logK logp(ci)=logK
又因为有
p ( x i ∣ μ , c i ) = ∏ k = 1 K p ( x i ∣ μ k ) c i k p(x_i|\mu,c_i) = \prod_{k=1}^Kp(x_i|\mu_k)^{c_{ik}} p(xiμ,ci)=k=1Kp(xiμk)cik

w h e r e               p ( x i ∣ u k ) = 1 2 π e − ( x i − u k ) 2 2 where\ \ \ \ \ \ \ \ \ \ \ \ \ p(x_i|u_k)=\frac {1}{\sqrt{2\pi}}e^{-\frac {{(x_i-u_k)}^2}{2}} where             p(xiuk)=2π 1e2(xiuk)2

上述公式是将给定所有的 μ \mu μ的条件下拆解成每个 μ k \mu_k μk乘积的形式, c i k c_{ik} cik的取值为0和1代表选择的意思

针对上述例子的第一个变量 x 1 x_1 x1此公式又可以写成
p ( x 1 ∣ μ , c 1 ) = p ( x 1 ∣ μ 1 ) c 11 p ( x 1 ∣ μ 2 ) c 12 p(x_1|\mu,c_1) = p(x_1|\mu_1)^{c_{11}}p(x_1|\mu_2)^{c_{12}} p(x1μ,c1)=p(x1μ1)c11p(x1μ2)c12
对于 E q − c i [ l o g p ( x i ∣ μ , c i ) ] E_{q-c_i}[logp(x_i|\mu,c_i)] Eqci[logp(xiμ,ci)]我们有
E q − c i [ l o g p ( x i ∣ μ , c i ) ] = ∑ k = 1 K c i k E q − c i [ p ( x i ∣ μ k ) ] E_{q-c_i}[logp(x_i|\mu,c_i)] = \sum_{k=1}^Kc_{ik}E_{q-c_i}[p(x_i|\mu_k)] Eqci[logp(xiμ,ci)]=k=1KcikEqci[p(xiμk)]
因为 E q − c i [ p ( x i ∣ μ k ) ] E_{q-c_i}[p(x_i|\mu_k)] Eqci[p(xiμk)]中只含有 u k u_k uk项因此除去 u k u_k uk项的其余项的期望为常数,得到以下的式子
E q − c i [ l o g p ( x i ∣ μ , c i ) ] = ∑ k = 1 K c i k E q − c i [ p ( x i ∣ μ k ) ] = ∑ k = 1 K c i k E q ( c 1 , c 2... c i − 1 , c i + 1 , . . . c n , u 1 , . . . u k − 1 , u k + 1 , . . . u K ) [ p ( x i ∣ μ k ) ] + ∑ k = 1 K c i k E q ( u k ; m k , s k 2 ) [ l o g p ( x i ∣ μ k ) ] E_{q-c_i}[logp(x_i|\mu,c_i)] =\sum_{k=1}^Kc_{ik}E_{q-c_i}[p(x_i|\mu_k)]=\sum_{k=1}^Kc_{ik}E_{q(c1,c2...c_{i-1},c_{i+1},...c_{n},u_1,...u_{k-1},u_{k+1},...u_K)}[p(x_i|\mu_k)]+ \sum_{k=1}^Kc_{ik}E_{q(u_k;m_k,s_k^2)}[logp(x_i|\mu_k)] Eqci[logp(xiμ,ci)]=k=1KcikEqci[p(xiμk)]=k=1KcikEq(c1,c2...ci1,ci+1,...cn,u1,...uk1,uk+1,...uK)[p(xiμk)]+k=1KcikEq(uk;mk,sk2)[logp(xiμk)]

= ∑ k = 1 K c i k E q ( u k ; m k , s k 2 ) [ − ( x i − μ k ) 2 2 ] + c o n s t =\sum_{k=1}^Kc_{ik}E_{q(u_k;m_k,s_k^2)}[-\frac{(x_i-\mu_k)^2}{2}]+const =k=1KcikEq(uk;mk,sk2)[2(xiμk)2]+const

= ∑ k = 1 K c i k [ E q ( u k ; m k , s k 2 ) [ μ k ] − E q ( u k ; m k , s k 2 ) [ μ k 2 ] ] + c o n s t =\sum_{k=1}^Kc_{ik}[E_{q(u_k;m_k,s_k^2)}[\mu_k]-E_{q(u_k;m_k,s_k^2)}[{\mu_k}^2]]+const =k=1Kcik[Eq(uk;mk,sk2)[μk]Eq(uk;mk,sk2)[μk2]]+const

对于变量 x 1 x_1 x1上述式子又可以写成
E q − c 1 [ l o g p ( x 1 ∣ μ , c 1 ) ] = c 11 [ E q ( u 1 ; m 1 , s 1 2 ) [ μ 1 ] − E q ( u 1 ; m 1 , s 1 2 ) [ μ 1 2 ] ] + c 12 [ E q ( u 2 ; m 2 , s 2 2 ) [ μ 2 ] − E q ( u 2 ; m 2 , s 2 2 ) [ μ 2 2 ] ] + c o n s t E_{q-c_1}[logp(x_1|\mu,c_1)] =c_{11}[E_{q(u_1;m_1,s_1^2)}[\mu_1]-E_{q(u_1;m_1,s_1^2)}[{\mu_1}^2]]+c_{12}[E_{q(u_2;m_2,s_2^2)}[\mu_2]-E_{q(u_2;m_2,s_2^2)}[{\mu_2}^2]]+const Eqc1[logp(x1μ,c1)]=c11[Eq(u1;m1,s12)[μ1]Eq(u1;m1,s12)[μ12]]+c12[Eq(u2;m2,s22)[μ2]Eq(u2;m2,s22)[μ22]]+const
又因为有 E q ( u k ; m k , s k 2 ) [ μ k ] = m k E_{q(u_k;m_k,s_k^2)}[\mu_k]=m_k Eq(uk;mk,sk2)[μk]=mk E q ( u k ; m k , s k 2 ) [ μ k 2 ] ] = v a r [ u k ] + [ E [ u k ] 2 ] = s k 2 + m k 2 E_{q(u_k;m_k,s_k^2)}[{\mu_k}^2]]=var[u_k]+[E[u_k]^2]=s_k^2+m_k^2 Eq(uk;mk,sk2)[μk2]]=var[uk]+[E[uk]2]=sk2+mk2因此我们得出最后的迭代公式
l o g q ∗ ( c i ; φ i ) ∝ ∑ k = 1 K c i k [ m k x i − ( s k 2 + m k 2 ) 2 ] logq^*(c_i;\varphi_i)\propto\sum_{k=1}^Kc_{ik}[\frac {m_kx_i-(s_k^2+m_k^2)}{2}] logq(ci;φi)k=1Kcik[2mkxi(sk2+mk2)]

q ∗ ( c i ; φ i ) ∝ e x p { ∑ k = 1 K c i k [ m k x i − ( s k 2 + m k 2 ) 2 ] } q^*(c_i;\varphi_i)\propto exp\{\sum_{k=1}^Kc_{ik}[\frac {m_kx_i-(s_k^2+m_k^2)}{2}]\} q(ci;φi)exp{k=1Kcik[2mkxi(sk2+mk2)]}

因为我们有 q ∗ ( c i = k ; φ i ) = φ i k q^*(c_i=k;\varphi_i)=\varphi_{ik} q(ci=k;φi)=φik所以有
q ∗ ( c i = k ; φ i ) = φ i k ∝ e x p { m k x i − ( s k 2 + m k 2 ) 2 } q^*(c_i=k;\varphi_i)=\varphi_{ik}\propto exp\{\frac {m_kx_i-(s_k^2+m_k^2)}{2}\} q(ci=k;φi)=φikexp{2mkxi(sk2+mk2)}
我们就得到了我们的更新公式,针对上述例子我们得到例子的更新公式如下
φ 11 = e x p { m 1 x 1 − ( s 1 2 + m 1 2 ) 2 } \varphi_{11} = exp\{\frac {m_1x_1-(s_1^2+m_1^2)}{2}\} φ11=exp{2m1x1(s12+m12)}

φ 12 = e x p { m 2 x 1 − ( s 2 2 + m 2 2 ) 2 } \varphi_{12} = exp\{\frac {m_2x_1-(s_2^2+m_2^2)}{2}\} φ12=exp{2m2x1(s22+m22)}

. . . . . . . ....... .......

φ 10 , 2 = e x p { m 2 x 10 − ( s 2 2 + m 2 2 ) 2 } \varphi_{10,2} = exp\{\frac {m_{2}x_{10}-(s_{2}^2+m_{2}^2)}{2}\} φ10,2=exp{2m2x10(s22+m22)}

接下来继续来求解 q ( μ k ; m k , s k ) q(\mu_k;m_k,s_k) q(μk;mk,sk)

同理,我们求解 q ∗ ( μ k ; m k , s k ) q^*(\mu_k;m_k,s_k) q(μk;mk,sk)
l o g q ∗ ( μ k ; m k , s k ) ∝ E q − u k [ l o g p ( c , u , x ) ] logq^*(\mu_k;m_k,s_k) \propto E_{q-u_k}[log p(c,u,x)] logq(μk;mk,sk)Equk[logp(c,u,x)]

E q − u k [ l o g p ( c , u , x ) ] = E q − u k [ l o g ( p ( u k ) p ( u − k ) p ( c ) p ( x ∣ μ , c ) ) ] E_{q-u_k}[log p(c,u,x)] = E_{q-u_k}[log(p(u_k)p(u_{-k})p(c)p(x|\mu,c))] Equk[logp(c,u,x)]=Equk[log(p(uk)p(uk)p(c)p(xμ,c))]

= E q − u k [ l o g p ( μ k ) + l o g p ( u − k ) + l o g p ( c ) + l o g p ( x ∣ μ , c ) ] =E_{q-u_k}[logp(\mu_k)+logp(u_{-k})+logp(c)+logp(x|\mu,c)] =Equk[logp(μk)+logp(uk)+logp(c)+logp(xμ,c)]

= E q − u k [ l o g p ( μ k ) ] + E q − u k [ l o g p ( u − k ) ] + E q − u k [ l o g p ( c ) ] + ∑ i = 1 n E q − u k [ l o g p ( x i ∣ μ , c i ) ] =E_{q-u_k}[logp(\mu_k)]+E_{q-u_k}[logp(u_{-k})]+E_{q-u_k}[logp(c)]+\sum_{i=1}^nE_{q-u_k}[logp(x_i|\mu,c_i)] =Equk[logp(μk)]+Equk[logp(uk)]+Equk[logp(c)]+i=1nEquk[logp(xiμ,ci)]

= l o g p ( u k ) + ∑ i = 1 n E q − u k [ l o g p ( x i ∣ μ , c i ) + c o n s t =logp(u_k)+\sum_{i=1}^nE_{q-u_k}[logp(x_i|\mu,c_i)+const =logp(uk)+i=1nEquk[logp(xiμ,ci)+const

针对上述例子中的 u 1 u_1 u1我们有如下表达式
E q − u 1 [ l o g p ( c , u , x ) ] = l o g p ( u 1 ) + E q − u 1 [ l o g p ( x 1 ∣ μ , c 1 ) + E q − u 1 [ l o g p ( x 2 ∣ μ , c 2 ) + . . . + E q − u 1 [ l o g p ( x 10 ∣ μ , c 10 ) + c o n s t E_{q-u_1}[log p(c,u,x)]=logp(u_1)+E_{q-u_1}[logp(x_1|\mu,c_1)+E_{q-u_1}[logp(x_2|\mu,c_2)+...+E_{q-u_1}[logp(x_{10}|\mu,c_{10})+const Equ1[logp(c,u,x)]=logp(u1)+Equ1[logp(x1μ,c1)+Equ1[logp(x2μ,c2)+...+Equ1[logp(x10μ,c10)+const
同样的我们将于 u k u_k uk无关的项都提取出来为常数项,接下来我们只需要求解以下式子即可
q ∗ ( μ k ; m k , s k ) ∝ e x p { l o g p ( u k ) + ∑ i = 1 n E q − u k [ l o g p ( x i ∣ μ , c i ) } q^*(\mu_k;m_k,s_k)\propto exp\{logp(u_k)+\sum_{i=1}^nE_{q-u_k}[logp(x_i|\mu,c_i)\} q(μk;mk,sk)exp{logp(uk)+i=1nEquk[logp(xiμ,ci)}
又因为有
l o g p ( u k ) = − μ k 2 2 σ 2 logp(u_k)=-\frac {\mu_k^2}{2\sigma^2} logp(uk)=2σ2μk2

∑ i = 1 n E q − u k [ l o g p ( x i ∣ μ , c i ) = ∑ i = 1 n E q − u k [ c i k l o g p ( x i ∣ μ k ) ] + c o n s t \sum_{i=1}^nE_{q-u_k}[logp(x_i|\mu,c_i)=\sum_{i=1}^nE_{q-u_k}[c_{ik}logp(x_i|\mu_k)]+const i=1nEquk[logp(xiμ,ci)=i=1nEquk[ciklogp(xiμk)]+const

= ∑ i = 1 n E q − u k [ c i k ] E q − u k [ l o g p ( x i ∣ μ k ) ] + c o n s t =\sum_{i=1}^nE_{q-u_k}[c_{ik}]E_{q-u_k}[logp(x_i|\mu_k)]+const =i=1nEquk[cik]Equk[logp(xiμk)]+const

对于上述式子 E q − u k [ c i k ] E_{q-u_k}[c_{ik}] Equk[cik]
E q − u k [ c i k ] = q ( c i = k ; φ i ) c i k = φ i k E_{q-u_k}[c_{ik}]=q(c_i=k;\varphi_i)c_{ik}=\varphi_{ik} Equk[cik]=q(ci=k;φi)cik=φik
对于上述例子中 u k u_k uk来说
E q − u 1 [ c i 1 ] = q ( c i = 1 ; φ i ) c i 1 = φ i 1 E_{q-u_1}[c_{i1}]=q(c_i=1;\varphi_i)c_{i1}=\varphi_{i1} Equ1[ci1]=q(ci=1;φi)ci1=φi1
继续上述式子的推导
∑ i = 1 n E q − u k [ l o g p ( x i ∣ μ , c i ) = ∑ i = 1 n E q − u k [ c i k ] E q − u k [ l o g p ( x i ∣ μ k ) ] + c o n s t \sum_{i=1}^nE_{q-u_k}[logp(x_i|\mu,c_i)=\sum_{i=1}^nE_{q-u_k}[c_{ik}]E_{q-u_k}[logp(x_i|\mu_k)]+const i=1nEquk[logp(xiμ,ci)=i=1nEquk[cik]Equk[logp(xiμk)]+const

= ∑ i = 1 n φ i k l o g p ( x i ∣ μ k ) + c o n s t =\sum_{i=1}^n\varphi_{ik}logp(x_i|\mu_k)+const =i=1nφiklogp(xiμk)+const

= ∑ i = 1 n φ i k [ − ( x i − μ k ) 2 2 ] + c o n s t =\sum_{i=1}^n\varphi_{ik}[-\frac{(x_i-\mu_k)^2}{2}]+const =i=1nφik[2(xiμk)2]+const

= ∑ i = 1 n φ i k x i μ k − ∑ i = 1 n φ i k μ k 2 2 + c o n s t =\sum_{i=1}^n\varphi_{ik}x_i\mu_k-\sum_{i=1}^n\varphi_{ik}\frac {{\mu_k}^2}{2}+const =i=1nφikxiμki=1nφik2μk2+const

因此有
q ∗ ( μ k ; m k , s k ) ∝ e x p { − μ k 2 2 σ 2 + ∑ i = 1 n φ i k x i μ k − 1 2 ∑ i = 1 n φ i k μ k 2 } q^*(\mu_k;m_k,s_k)\propto exp\{-\frac {\mu_k^2}{2\sigma^2}+\sum_{i=1}^n\varphi_{ik}x_i\mu_k-\frac {1}{2}\sum_{i=1}^n\varphi_{ik}{\mu_k}^2\} q(μk;mk,sk)exp{2σ2μk2+i=1nφikxiμk21i=1nφikμk2}

∝ e x p { ∑ i = 1 n ( φ i k x i ) μ k − 1 2 ( σ 2 + ∑ i = 1 n φ i k ) μ k 2 } \propto exp\{\sum_{i=1}^n(\varphi_{ik}x_i)\mu_k-\frac {1}{2}(\sigma^2+\sum_{i=1}^n\varphi_{ik}){\mu_k}^2\} exp{i=1n(φikxi)μk21(σ2+i=1nφik)μk2}

因为 q ( u k ) q(u_k) q(uk)是高斯分布,我们按照高斯分布的形式 1 2 π e x p ( x − ϵ ) 2 2 s k 2 \frac {1}{\sqrt{2\pi}}exp^{\frac {(x-\epsilon)^2}{2s_k^2}} 2π 1exp2sk2(xϵ)2我们对指数部分进行拆解,
− ( x − m k ) 2 2 s k 2 = − ( x 2 − 2 x m k + m k 2 ) 2 s k 2 {\frac {-(x-m_k)^2}{2s_k^2}}=\frac {-(x^2-2xm_k+m_k^2)}{2s_k^2} 2sk2(xmk)2=2sk2(x22xmk+mk2)

= − x 2 2 s k 2 + 2 x m k 2 s k 2 − m k 2 2 s k 2 =-\frac {x^2}{2s_k^2} +\frac{2xm_k}{2s_k^2}-\frac{m_k^2}{2s_k^2} =2sk2x2+2sk22xmk2sk2mk2

− 1 2 ( σ 2 + ∑ i = 1 n φ i k ) = − 1 2 s k 2 -\frac {1}{2}(\sigma^2+\sum_{i=1}^n\varphi_{ik})=-\frac {1}{2s_k^2} 21(σ2+i=1nφik)=2sk21

求得
s k = 1 σ 2 + ∑ i = 1 n φ i k s_k=\frac {1}{{σ^2}+∑_{i=1}^nφ_{ik}} sk=σ2+i=1nφik1

对于 m k m_k mk来说
m k s k = ∑ i = 1 n ( φ i k x i ) \frac {m_k}{s_k} =\sum_{i=1}^n(\varphi_{ik}x_i) skmk=i=1n(φikxi)

m k = ∑ i = 1 n φ i k x i σ 2 + ∑ i n φ i k m_k =\frac {∑_{i=1}^n φ_{ik}x_i}{{σ^2}+∑_i^nφ_{ik}} mk=σ2+inφiki=1nφikxi

例如我的第一次对第1个分布的均值,方差的迭代就应该如下所示:
s 1 = 1 100 + φ 1 , 1 + φ 2 , 1 + φ 3 , 1 + . . . . , + φ 10 , 1 s_1 = \frac {1}{{100}+\varphi_{1,1}+\varphi_{2,1}+\varphi_{3,1}+....,+\varphi_{10,1}} s1=100+φ1,1+φ2,1+φ3,1+....,+φ10,11

m 1 = φ 1 , 1 x 1 + φ 2 , 1 x 2 + . . . , + φ 10 , 1 x 10 100 + φ 1 , 1 + φ 2 , 1 + φ 3 , 1 + . . . . , + φ 10 , 1 m_1 = \frac {\varphi_{1,1}x_{1}+\varphi_{2,1}x_{2}+...,+\varphi_{10,1}x_{10}}{{100}+\varphi_{1,1}+\varphi_{2,1}+\varphi_{3,1}+....,+\varphi_{10,1}} m1=100+φ1,1+φ2,1+φ3,1+....,+φ10,1φ1,1x1+φ2,1x2+...,+φ10,1x10

最后我们得到了所有参数的更新公式,总结求解过程为:

  • 随机初始化所有参数 φ , m , s \mathbf \varphi, \mathbf m, \mathbf s φ,m,s
  • 更新每一个参数参数$\varphi_i $
  • 更新每一个参数 s k s_k sk
  • 更新每一个参数 m k m_k mk按照(4)式计算ELBO, 如果ELBO收敛则结束,否则返回第1步
代码实现:

根据上述算法的描述,可以写出实现代码,但第5步可以不需要计算ELBO, 直接迭代n步之后结束即可,具体代码如下:

def solve(x, k, sigma, epoch=40):
    """
    x: 输入数据
    k: 超参数k,c_i的维度,在业务CASE中等于用户数
    sigma: 超参数,需要人工调整
    """
    n = len(x)
    phis = np.random.random([n, k])
    mk = np.random.random([k])
    sk = np.random.random([k])
    for _ in range(epoch):
        for i in range(n):
            phi_i_k = []
            for _k in range(k):
                #根据公式(6)更新参数phi_ik
                phi_i_k.append(np.exp(mk[_k]*x[i] - (sk[_k]**2 + mk[_k]**2)/2))
            sum_phi = sum(phi_i_k)
            phi_i_k = [phi/sum_phi for phi in phi_i_k]
            phis[i] = phi_i_k
        den = np.sum(phis, axis=0) + 1/(sigma**2)
        #根据公式(10)更新m_k
        mk = np.matmul(x, phis)/den
        #根据公式(11)更新s_k
        sk = np.sqrt(1/den)
    return mk, sk, phis

mk, sk, phis = solve(x,k,sigma)

输入第二节中生成的数据x和超参数后,求解得到m:[17.08924954 4.86667515], 对比第二节的真实参数u十分接近。从图像上, 根据求解出来的参数 m , s , φ \mathbf m, \mathbf s, \varphi m,s,φ,模拟采样数据,得到的数据分布也与真实数据分布十分接近(黄色为真实数据,蓝色模拟采样数据)

在这里插入图片描述

完整代码
import numpy as np
import random
import matplotlib.pyplot as plt
import seaborn as sns

sigma = 10
k = 2
n = 10

def gen_data(sigma, k, n):
    #获取u_k
    u_k = np.random.normal(0, sigma, k)
    x = []
    c = []
    for i in range(n):
        ci = random.choice(range(k))
        c.append(ci)
        u = u_k[ci]
        x.append(np.random.normal(u, 1))
    return x, u_k, c

x, u_k, c = gen_data(sigma,k,n)

print("x:"+str(x))
print("u_k:"+str(u_k))
print("c:"+str(c))

sns.distplot(x, hist=False,color='y')
sns.distplot(x,color='y')
plt.show()

print('**'*100)

def solve(x, k, sigma, epoch=40):
    """
    x: 输入数据
    k: 超参数k,c_i的维度,在业务CASE中等于用户数
    sigma: 超参数,需要人工调整
    """
    n = len(x)
    phis = np.random.random([n, k])
    mk = np.random.random([k])
    sk = np.random.random([k])
    for _ in range(epoch):
        for i in range(n):
            phi_i_k = []
            for _k in range(k):
                #根据公式(6)更新参数phi_ik
                phi_i_k.append(np.exp(mk[_k]*x[i] - (sk[_k]**2 + mk[_k]**2)/2))
            sum_phi = sum(phi_i_k)
            phi_i_k = [phi/sum_phi for phi in phi_i_k]
            phis[i] = phi_i_k
        den = np.sum(phis, axis=0) + 1/(sigma**2)
        #根据公式(10)更新m_k
        mk = np.matmul(x, phis)/den
        #根据公式(11)更新s_k
        sk = np.sqrt(1/den)
    return mk, sk, phis

mk, sk, phis = solve(x,k,sigma)
n = 10 # number of sample to be drawn
samples = []
for i in range(n): # iteratively draw samples
    Z = np.random.choice([0,1]) # latent variable
    samples.append(np.random.normal(mk[Z], sk[Z], 1))

sns.distplot(x, hist=False,color='y')
sns.distplot(x,color='y')
sns.distplot(samples, hist=False,color='b')
sns.distplot(samples,color='b')
plt.show()
print("mk:"+str(mk))
print("sk:"+str(sk))
print("phis:"+str(phis))
pred = []
for p in phis:
    pred.append(np.argmax(p))
print(pred)

上个模拟数据的结果:

m m m:[17.08924954 4.86667515]
s s s:[0.40790851 0.49937617]

c l a s s class class:[1, 0, 1, 0, 1, 0, 1, 0, 0, 0]

拓展:如何使用更复杂的后验分布簇来降低VI方法的bias

变分分布是用来替代真实后验分布的,两者的差异越大,后验推断的系统偏差就会越大。有研究结果表明,变分后验分布簇的选择对变分推断效果的影响非常大。

经典的 VI,会基于简单的平均场(mean-fifiled)假设,用可分解的高斯分布或者一些简单结构的分布来作为变分分布;现在的 VI,需要解决的是数据规模更大、维度更高的问题,经典 VI 的变分分布难以满足。因此,最近几年有一系列工作来研究如何构造一系列更加复杂且方便计算的复杂后验分布来解决这一问题。

在这里插入图片描述

Copula方法

大多数的 VI 方法都基于 Mean-Field 的思路,假设变分后验分布中隐变量之间相互独立,这个假设太强,对结果有一定的影响。

NIPS 2015 一篇 David M. Blei 组的工作 Copula Variational Inference 尝试用统计学的经典方法 Copula 来解决 MF 中隐变量的独立假设问题。这篇工作的动机非常简单,就是找到一种既考虑隐变量之间的关联性同时也容易进行大规模计算的方法。思路如下:
q ( z ; θ , η ) = [ ∏ i = 1 d q ( z i ; θ ) c ( Q ( z i ; θ ) , . . . . , Q ( z d ; θ ) ; η ) ] q(z;\theta,\eta)=[\prod_{i=1}^dq(z_i;\theta)c(Q(z_i;\theta),....,Q(z_d;\theta);\eta)] q(z;θ,η)=[i=1dq(zi;θ)c(Q(zi;θ),....,Q(zd;θ);η)]
其中公式中的前半部分是 Mean-Field,而后半部分正是所谓的 Copula。

辅助变量法

在这里插入图片描述

辅助变量法的思路比较简单,它认为隐变量 z 背后还有隐变量 w,是一种层次化建模的思想。即:
q ( z ) = ∫ q ( z ∣ w ) q ( w ) d w q(z)=\int q(z|w)q(w)dw q(z)=q(zw)q(w)dw

评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

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

余额充值