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

k1、变分推断的算法

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

X：观测变量，Z:隐变量+参数

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

$$p(Z|X)=\frac{p(X,Z)}{{\int }_{Z}p(X|Z)p(Z)dZ}$$

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

进行简单的变换我们可以得到
$$p(X)=\frac{p(X,Z)}{p(Z|X)}$$
两边同时取对数可以得到：
$$logP(X)=log\frac{p(X,Z)}{p(Z|X)}$$

$$=logp(X,Z)-logp(Z|X)$$

$$=log\frac{p(X,Z)}{q(Z)}-log\frac{p(Z|X)}{q(Z)}$$

$$=logp(X,Z)-logq(Z)-log\frac{p(Z|X)}{q(Z)}$$

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

1.1公式化简

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

左边：

对于左边的式子$logp(X)$求q(Z)的期望得：
$$E(logp(X))={\int }_{Z}logp(X)q(Z)dZ=logp(X){\int }_{Z}q(Z)dZ$$
由于${\int }_{Z}q(Z)dZ=1$故最终的$E(logp(x))=logp(X)$

右边：

$$logp(X,Z)-logq(Z)-log\frac{p(Z|X)}{q(Z)}=\underset{ELBO}{\underbrace{{\int }_{Z}logp(X,Z)q(Z)dZ-{\int }_{Z}logq(Z)q(Z)dZ}}+\underset{KL}{\underbrace{{\int }_Z-log\frac{p(Z|X)}{q(Z)}q(Z)dZ}}$$

$$=\underset{ELBO}{\underbrace{{\int }_{Z}log\frac{p(X,Z)}{q(Z)}q(Z)dZ}}+\underset{KL}{\underbrace{{\int }_Z-log\frac{p(Z|X)}{q(Z)}q(Z)dZ}}$$

最终化简的形式如下：
$$logp(X）=\underset{ELBO}{\underbrace{{\int }_{Z}log\frac{p(X,Z)}{q(Z)}q(Z)dZ}}+\underset{KL}{\underbrace{{\int }_Z-log\frac{p(Z|X)}{q(Z)}q(Z)dZ}}$$

其中$ELBO$也是被记为$L(q)$

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

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

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

我们想要求得一个 $q(Z)\approx p(Z|X)$。也就是
$$q(Z)=argma{x}_{q(z)}L(q)$$

1.2模型求解（经典平均场理论）

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

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

通过平均场理论我们就把原先高维的隐变量转化成低维的独立的隐变量相乘，极大简化的积分的复杂度
$$L(q)={\int }_{Z}logp(X,Z)q(Z)dZ-{\int }_{Z}logq(Z)q(Z)dZ$$

$$={\int }_{Z}logp(X,Z)\prod _{i=1}^M{q}_i(z_i)dZ-{\int }_Z\sum _{i=1}^{M}logq(Z)\prod _{i=1}^M{q}_i(z_i)dZ$$

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

part 1：

$${\int }_{Z}logp(X,Z)q(Z)dZ={\int }_{Z}logp(X,Z)\prod _{i=1}^M{q}_i(z_i)dZ$$

我们将整个集合的Z拆解成一个个z进行积分
$$={\int }_{z_1}{\int }_{z_2}...{\int }_{z_M}{q}_i(z_i)log(p(X,Z))d{z}_{1}d{z}_2...d{z}_M$$
关键的地方来了，我们这里除了第$Z_j$个隐变量讲其他项的都积分掉，这样就可以留下含有$Z_j$项的式子，这样我们后面就可以通过CAVI(坐标上升法)来对每个参数进行迭代更新，
$$={\int }_{z_j}{q}_j(z_j)[{\int }_{z_1}{\int }_{z_2}...{\int }_{z_M}{q}_i(z_i)log(p(X,Z))d{z}_{1}d{z}_2...d{z}_M]d_j$$

$$={\int }_{z_j}{q}_j(z_j)E_{{\prod }_{i\ne j}^M{q}_i(z_i)}[logp(X,Z)]d{z}_j$$

part 2:

$${\int }_Z\sum _{i=1}^{M}logq(Z)\prod _{i=1}^M{q}_i(z_i)dZ={\int }_Z\prod _{i=1}^M{q}_i(z_i)dZ[log{q}_1(z_1)+log{q}_2(z_2)+...+log{q}_M(z_M)]dZ$$

我们如何计算呢？现在我们通过简单的式子来化简出通向公式，现在设定$M=2$有
$${\int }_{z_1}{\int }_{z_2}[log{q}_1(z_1)+log{q}_2(z_2)]q_1(z_1)q_2(z_2)d{z}_{1}d{z}_2={\int }_{z_1}{\int }_{z_2}{q}_1(z_1)q_2(z_2)log{q}_1(z_1)d{z}_{1}d{z}_2+{\int }_{z_1}{\int }_{z_2}{q}_1(z_1)q_2(z_2)log{q}_2(z_2)d{z}_{1}d{z}_2$$

$$={\int }_{z_1}{q}_1(z_1)log{q}_1(z_1)\underset{1}{\underbrace{{\int }_{z_2}{q}_2(z_2)d{z}_2}}d{z}_1+{\int }_{z_2}{q}_2(z_2)log{q}_2(z_2)\underset{1}{\underbrace{{\int }_{z_1}{q}_1(z_1)d{z}_1}}d{z}_2$$

$$=\sum _{i=1}^2{\int }_{z_i}{q}_i(z_i)log{q}_i(z_i)d{z}_i$$

故可以求得通向公式
$${\int }_Z\sum _{i=1}^{M}logq(Z)\prod _{i=1}^M{q}_i(z_i)dZ=\sum _{i=1}^M{\int }_{z_i}{q}_i(z_i)log{q}_i(z_i)d{z}_i$$
因为我们仅仅只关注第 j 项因为第j项是我们需要进行迭代的一项，所以除了j项其余项都为常数
$$part2={\int }_{z_j}{q}_i(z_i)log{q}_i(z_i)d{z}_j+C$$

为了近一步计算我们将上述$part1$式子$E_{{\prod }_{i\ne j}^M{q}_i(z_i)}[logp(X,Z)]$表示为
$$E_{{\prod }_{i\ne j}^M{q}_i(z_i)}[logp(X,Z)]=log\hat{p}(X,z_j)$$
为什么要这样做呢，目的是为了可以构造出$L(q)=-KL$后面的式子，求出我们每个$Z_j$分布需要逼近的真实的分布，我们可以对上述式子求逆得到伪真实分布$\hat{p}(X,z_j)$为什么可以叫做伪真实分布因为它们之间只是相差了一个常数项
$$\hat{p}(X,z_j)=ex{p}^{E_{{\prod }_{i\ne j}^M{q}_i(z_i)}[logp(X,Z)]d{z}_j}$$
那么$part1$的式子就可以写做
$$={\int }_{z_j}{q}_j(z_j)log\hat{p}(X,z_j)d{z}_j$$
这样我们就可以求得最终的$L(q)$的表达式
$$L(q_j)=part1-part2={\int }_{z_j}{q}_j(z_j)log\hat{p}(X,z_j)d{z}_j-{\int }_{z_j}{q}_i(z_j)log{q}_i(z_j)d{z}_j+C$$

$$={\int }_{z_j}{q}_j(z_j)log\frac{\hat{p}(X,z_j)}{q_i(z_j)}d{z}_j+C=-KL$$

$$={\int }_{z_j}{q}_j(z_j)log\hat{p}(X,z_j)d{z}_j+C$$

因此最后我们要求解$argma{x}_{q(z)}L(q_j)$就是求解等价于$argmi{n}_{q(z)}KL(q_j||\hat{p}(x,z_j))$ ，有kl散度的定义我们可以知道只有当$q_j$等于$\hat{p}(x,z_j)$时可以达到最优解，对此我们就可以得到每个隐变量的迭代更新的式子
$$q_i(z_i)=ex{p}^{E_{{\prod }_{i\ne j}^M{q}_i(z_i)}[log\hat{p}(X,Z)]d{z}_j}$$

$$log{q}_i(z_i)=E_{{\prod }_{i\ne j}^M{q}_i(z_i)}[log\hat{p}(X,Z)]d{z}_j$$

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

1.3一个业务的例子

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

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

• 某个用户$c_k$登录系统
• 用户产生一条使用时长增量记录$x_i$

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

${\mu }_k$代表的是我的观测数据$x_i$所属的分布的均值
$${\mu }_k\sim N(0,{\sigma }^2)k=1,....K$$
我们可以写出$u_k$的概率密度函数
$$p(u_k)=\frac{1}{\sqrt{2\pi {\sigma }^2}}{e}^{-\frac{u_k^2}{2{\sigma }^2}}$$
$c_i$服从类别分布
$$c_i\sim Categorical(\frac{1}{k},\frac{1}{k},...\frac{1}{k})i=1,.....,n$$
$p(c_i)$的概率密度函数如下：
$$p(c_i)=\frac{1}{k}$$
$c_i$代表的是第i条观测数据$x_i$是来自于哪个分布,例如：当K=6时$ci=[0,0,1,0,0,0]$的向量,代表第i条数据$x_i$属于第3个分布
$$x_i|c_i\sim N(c_i,\mu ,1)i=1,....,n$$

值得注意的是由于我们的假设$p({\mu }_k)$和$p(c_i)$是相互独立的先验概率，后面我们需要通过一个我们定义的$q(z)$分布去拟合这些真实的分布$p(z|x)$注意这里是给定x条件的z的分布，这样的目的是我们可以把原先的先验分布转化成后验分布，那么去拟合后验分布有什么好处呢，就举一个简单的例子，对于判断给定的变量x是来自哪个分布的$u_k$如果按照先验分布去拟合，最后得到的结果就是对于每个变量x来自于每个$u_k$的概率都是一样的，这显然是毫无意义的，我们的目的就是为了让模型可以帮助我们把变量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$:[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\in k}$的距离，选择距离最近的$ i=argmin_{i} dis(u_{i \in k}, x’)$赋予该条数据作为其用户id即可。

接下来我们来求我们需要求的隐变量$u_k,c_i$所服从的分布对应的参数的值，这样当我们求解得出这些分布的参数后，我们就可以通过这些分布参数来求解这些隐变量最有可能的取值。从而推测出这些隐变量的值，对于上述隐变量我们需要求[${\phi }_{1,0},...,{\phi }_{10,0},{\phi }_{1,1},...,{\phi }_{10,1},m_0,s_0,m_1,s_1$]总共24个参数

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

${\phi }_{ik}$代表的是属于该分布的概率，对应的因为我的c是满足分类分布$c_{ik}$的取值范围只能是0或1，$，q(c_i=2;{\phi }_i)={\phi }_{i2}，q(c_i=1;{\phi }_i)={\phi }_{i1}$

先求解$q(c_i;{\phi }_i)$

有上述证明的$log{q}_j^{\ast }(z_j)\propto E_{{\prod }_{i\ne j}^M{q}_i(z_i)}[logp(z,x)]$(只有通过使用平均随机场时才可以用这个公式,

这里我们回到本例例子，求$q^{\ast }(c_i;{\phi }_i)$,带入公式得
$$log{q}^{\ast }(c_i;{\phi }_i)\propto E_{{\prod }_{j\ne i}^M{q}_j(z_j)}[logp(c,u,x)]$$
这里我们的$q_j(z_j)$代表$q_1(c_1),q_2(c_2),....q_{10}(c_{10}),q_{11}(u_1)q_{12}(u_2)$共计12个z
$$E_{{\prod }_{j\ne i}^M{q}_j(z_j)}[logp(c,u,x)]=E_{q-c_i}[log[p(\mu )p(c)p(x|\mu ,c)]]$$

$$=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})]]$$

$$=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})]$$

$$=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})]$$

从上述推导中，仅仅只有$E_{q-c_i}[logp(x_i|\mu ,c_i)]$和$logp(c_i)$是依赖于$c_i$的项，其余项都可以提出为常数项，因此我们就可以更新我们的正比公式
$$log{q}^{\ast }(c_i;{\phi }_i)\propto E_{q-c_i}[logp(x_i|\mu ,c_i)]+logp(c_i)$$
通过先验，我们可以得到
$$logp(c_i)=-logK$$
又因为有
$$p(x_i|\mu ,c_i)=\prod _{k=1}^{K}p(x_i|{\mu }_k{)}^{c_{ik}}$$

$$wherep(x_i|u_k)=\frac{1}{\sqrt{2\pi }}{e}^{-\frac{{(x_i-u_k)}^2}{2}}$$

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

针对上述例子的第一个变量$x_1$此公式又可以写成
$$p(x_1|\mu ,c_1)=p(x_1|{\mu }_1{)}^{c_{11}}p(x_1|{\mu }_2{)}^{c_{12}}$$
对于$E_{q-c_i}[logp(x_i|\mu ,c_i)]$我们有
$$E_{q-c_i}[logp(x_i|\mu ,c_i)]=\sum _{k=1}^K{c}_{ik}{E}_{q-c_i}[p(x_i|{\mu }_k)]$$
因为$E_{q-c_i}[p(x_i|{\mu }_k)]$中只含有$u_k$项因此除去$u_k$项的其余项的期望为常数，得到以下的式子
$$E_{q-c_i}[logp(x_i|\mu ,c_i)]=\sum _{k=1}^K{c}_{ik}{E}_{q-c_i}[p(x_i|{\mu }_k)]=\sum _{k=1}^K{c}_{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}^K{c}_{ik}{E}_{q(u_k;m_k,s_k^2)}[logp(x_i|{\mu }_k)]$$

$$=\sum _{k=1}^K{c}_{ik}{E}_{q(u_k;m_k,s_k^2)}[-\frac{(x_i-{\mu }_k{)}^2}{2}]+const$$

$$=\sum _{k=1}^K{c}_{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$$

对于变量$x_1$上述式子又可以写成
$$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$$
又因为有$E_{q(u_k;m_k,s_k^2)}[{\mu }_k]=m_k$和$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$因此我们得出最后的迭代公式
$$log{q}^{\ast }(c_i;{\phi }_i)\propto \sum _{k=1}^K{c}_{ik}[\frac{m_k{x}_i-(s_k^2+m_k^2)}{2}]$$

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=1∑K​cik​[2mk​xi​−(sk2​+mk2​)​]}

因为我们有$q^{\ast }(c_i=k;{\phi }_i)={\phi }_{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​)=φik​∝exp{2mk​xi​−(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{2m1​x1​−(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{2m2​x1​−(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{2m2​x10​−(s22​+m22​)​}

接下来继续来求解$q({\mu }_k;m_k,s_k)$

同理，我们求解$q^{\ast }({\mu }_k;m_k,s_k)$
$$log{q}^{\ast }({\mu }_k;m_k,s_k)\propto E_{q-u_k}[logp(c,u,x)]$$

$$E_{q-u_k}[logp(c,u,x)]=E_{q-u_k}[log(p(u_k)p(u_{-k})p(c)p(x|\mu ,c))]$$

$$=E_{q-u_k}[logp({\mu }_k)+logp(u_{-k})+logp(c)+logp(x|\mu ,c)]$$

$$=E_{q-u_k}[logp({\mu }_k)]+E_{q-u_k}[logp(u_{-k})]+E_{q-u_k}[logp(c)]+\sum _{i=1}^n{E}_{q-u_k}[logp(x_i|\mu ,c_i)]$$

$$=logp(u_k)+\sum _{i=1}^n{E}_{q-u_k}[logp(x_i|\mu ,c_i)+const$$

针对上述例子中的$u_1$我们有如下表达式
$$E_{q-u_1}[logp(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$$
同样的我们将于$u_k$无关的项都提取出来为常数项，接下来我们只需要求解以下式子即可
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=1∑n​Eq−uk​​[logp(xi​∣μ,ci​)}
又因为有
$$logp(u_k)=-\frac{{\mu }_k^2}{2{\sigma }^2}$$

$$\sum _{i=1}^n{E}_{q-u_k}[logp(x_i|\mu ,c_i)=\sum _{i=1}^n{E}_{q-u_k}[c_{ik}logp(x_i|{\mu }_k)]+const$$

$$=\sum _{i=1}^n{E}_{q-u_k}[c_{ik}]E_{q-u_k}[logp(x_i|{\mu }_k)]+const$$

对于上述式子$E_{q-u_k}[c_{ik}]$有
$$E_{q-u_k}[c_{ik}]=q(c_i=k;{\phi }_i)c_{ik}={\phi }_{ik}$$
对于上述例子中$u_k$来说
$$E_{q-u_1}[c_{i1}]=q(c_i=1;{\phi }_i)c_{i1}={\phi }_{i1}$$
继续上述式子的推导
$$\sum _{i=1}^n{E}_{q-u_k}[logp(x_i|\mu ,c_i)=\sum _{i=1}^n{E}_{q-u_k}[c_{ik}]E_{q-u_k}[logp(x_i|{\mu }_k)]+const$$

$$=\sum _{i=1}^n{\phi }_{ik}logp(x_i|{\mu }_k)+const$$

$$=\sum _{i=1}^n{\phi }_{ik}[-\frac{(x_i-{\mu }_k{)}^2}{2}]+const$$

$$=\sum _{i=1}^n{\phi }_{ik}{x}_i{\mu }_k-\sum _{i=1}^n{\phi }_{ik}\frac{{{\mu }_k}^2}{2}+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=1∑n​φik​xi​μk​−21​i=1∑n​φik​μk​2}

∝ 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=1∑n​(φik​xi​)μk​−21​(σ2+i=1∑n​φik​)μk​2}

因为$q(u_k)$是高斯分布，我们按照高斯分布的形式$\frac{1}{\sqrt{2\pi }}ex{p}^{\frac{(x-ϵ{)}^2}{2s_k^2}}$我们对指数部分进行拆解，
$$\frac{-(x-m_k{)}^2}{2s_k^2}=\frac{-(x^2-2x{m}_k+m_k^2)}{2s_k^2}$$

$$=-\frac{x^2}{2s_k^2}+\frac{2x{m}_k}{2s_k^2}-\frac{m_k^2}{2s_k^2}$$

$$-\frac{1}{2}({\sigma }^2+\sum _{i=1}^n{\phi }_{ik})=-\frac{1}{2s_k^2}$$

求得
$$s_k=\frac{1}{{\sigma }^2+{\sum }_{i=1}^n{\phi }_{ik}}$$

对于$m_k$来说
$$\frac{m_k}{s_k}=\sum _{i=1}^n({\phi }_{ik}{x}_i)$$

$$m_k=\frac{{\sum }_{i=1}^n{\phi }_{ik}{x}_i}{{\sigma }^2+{\sum }_i^n{\phi }_{ik}}$$

例如我的第一次对第1个分布的均值，方差的迭代就应该如下所示：
$$s_1=\frac{1}{100+{\phi }_{1,1}+{\phi }_{2,1}+{\phi }_{3,1}+....,+{\phi }_{10,1}}$$

$$m_1=\frac{{\phi }_{1,1}{x}_1+{\phi }_{2,1}{x}_2+...,+{\phi }_{10,1}{x}_{10}}{100+{\phi }_{1,1}+{\phi }_{2,1}+{\phi }_{3,1}+....,+{\phi }_{10,1}}$$

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

• 随机初始化所有参数$\phi ,\mathbf{m},\mathbf{s}$
• 更新每一个参数参数$\varphi_i $
• 更新每一个参数$s_k$
• 更新每一个参数$m_k$按照(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十分接近。从图像上， 根据求解出来的参数 $\mathbf{m},\mathbf{s},\phi$,模拟采样数据，得到的数据分布也与真实数据分布十分接近(黄色为真实数据，蓝色模拟采样数据)

完整代码

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$:[17.08924954 4.86667515]
$s$:[0.40790851 0.49937617]

$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;\theta ,\eta )=[\prod _{i=1}^{d}q(z_i;\theta )c(Q(z_i;\theta ),....,Q(z_d;\theta );\eta )]$$
其中公式中的前半部分是 Mean-Field，而后半部分正是所谓的 Copula。

辅助变量法

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