前言

参考：

• 《统计学习》—李航（蓝皮）

期望最大化算法（EM算法）

引入

简介

概率模型有时既含有观测变量（Observable Variable）又含有隐变量或潜在变量（Latent Variable）

e.g.购物网站，用户浏览商品行为包括点击促销广告/浏览高价商品/购买了一定频次，这些行为的关键点都在于平台只能看到这些行为数据,但不知道真实的用户类型，那么前面的这些数据便为观测变量，而真实的用户类型（用户画像）便为隐变量

EM算法就是含有隐变量的概率模型的参数极大似然估计法，或极大后验概率估计法

三硬币模型

介绍

假设有三枚硬币，分别记作A，B，C。这些硬币正面出现的概率分别是\pi,p,q，进行如下的掷硬币试验：
先掷硬币A，根据其结果选出硬币B或C
之后掷选出的硬币。结果正面记1反面记0
独立重复n次（这里n=10）
观测结果如下：
1，1，0，1，0，0，1，0，1，1
假设只能观测到掷硬币的结果，不能观测掷硬币的过程，问如何估计三硬币正面出现的概率，即三硬币模型的参数

解：
三硬币模型可写为
::: align-center
\begin{aligned} &P(y|\theta)=\sum_{z}P(y,z|\theta)\\ &=\sum_{z}P(z|\theta)P(y|z,\theta)\\ &=\pi p^y(1-p)^{1-y}+(1-\pi)q^y(1-q)^{1-y} \end{aligned}

:::
这里随机变量Y是观测变量，表示一次观测的结果1或0，随机变量Z是隐变量，Y和Z连在一起称为完全数据（Complete-Data），观测数据Y又称为不完全数据（Incomplete-Data），\theta=(\pi,p,q)是模型参数

• P(y,z|\theta)：在参数条件下观察到y,z的联合概率
• P(z|\theta)：在参数条件下z出现的概率
• P(y|z,\theta)：在参数和z出现条件下y出现的概率

将观测数据表示为Y=(Y_1,Y_2,\cdots,Y_n)^T，未观测数据表示为Z=(Z_1,Z_2,\cdots,Z_n)^T，则观测数据的似然函数为
::: align-center
P(Y|\theta)=\prod_{j=1}^n[\pi p^{y_j}(1-p)^{1-y_j}+(1-\pi)q^{y_j}(1-q)^{1-y_j}]

:::
考虑求解模型参数
::: align-center
\hat{\theta}=\arg\min_{\theta}\log P(Y|\theta)

:::
这个问题没有解析解，若强行求导数求解，各个变量间会相互耦合导致组合爆炸，这会使导数计算异常困难，所以要考虑使用其他方法求解

三硬币问题的EM算法

下面给出针对以上问题的EM算法：
EM算法首先选取参数的初值，记作\theta^{(0)}=(\pi^{(0)},p^{(0)},q^{(0)})，然后通过下面的步骤迭代计算参数的估计值，直至收敛为止——假设第i次迭代参数的估计值为\theta^{(i)}=(\pi^{(i)},p^{(i)},q^{(i)})，EM算法的第i+1次迭代如下：

E步（期望步）：计算在模型参数\theta^{(i)}下观测数据y_j来自掷硬币B的概率（即计算在当前参数下属于各隐变量的概率）
::: align-center
\mu_{j}^{(i+1)}=\frac{\pi p^{y_j}(1-p)^{1-y_j}}{\pi p^{y_j}(1-p)^{1-y_j}+(1-\pi)q^{y_j}(1-q)^{1-y_j}}

:::
这里其实使用到了贝叶斯公式，我们在计算P(Z=1|Y,\theta^{(i)})

之后计算在计算隐变量Z在条件分布P(Y,Z|\theta)下的期望，即
::: align-center
\begin{aligned} &Q(\theta,\theta^{(i)})=E_Z[\log(Y,Z|\theta)|Y,\theta^{(i)}]\\ &=\sum_Z\log P(Y,Z|\theta)P(Z|Y,\theta^{(i)}) \end{aligned}

:::

这里对于P(Y,Z|\theta)，我们考虑抛硬币的结果Z=0,1（伯努利0-1分布），那么我们就可以给出P(Y,Z|\theta)的似然函数为
::: align-center
\prod_{j=1}^n\{[\pi p^{y_j}(1-p)^{1-y_j}]^{z_j}[(1-\pi)q^{y_j}(1-q)^{1-y_j}]^{1-z_j}\}

:::
注：由P(Y,Z|\theta)=P(Y|Z,\theta)P(Z|\theta)和Z=0,1，分拆为两种情况，写为分段函数后合并结果（将z_j的分情况讨论转化为上面的幂次项）写出

取对数后，对于单次试验j，它的对数似然函数为
::: align-center
\begin{aligned} &\log P(y_j,z_j|\theta)\\ &=z_j[\log\pi+y_j\log p+(1-y_j)\log (1-p)]\\ &+(1-z_j)[\log (1-\pi)+y_j\log q+(1-y_j)\log (1-q)] \end{aligned}

:::
然后考虑z_j的条件期望，它实际上就是\mu_j，因为z_j只可能为0或1，这里的\mu_j是z_j=1的条件概率，故立即得结果

::: align-center
\begin{aligned} &Q(\theta,\theta^{(i)})\\ &=\sum_{j=1}^n\mu_j[\log\pi+y_j\log p+(1-y_j)\log (1-p)]\\ &+(1-\mu_j)[\log (1-\pi)+y_j\log q+(1-y_j)\log (1-q)]\\ \end{aligned}

:::

M步（最大化步）：计算模型参数的新估计值
::: align-center
\begin{aligned} &\pi^{(i+1)}=\frac{1}{n}\sum_{j=1}^n\mu_j^{(i+1)}\\ &p^{(i+1)}=\frac{\sum_{j=1}^n\mu_j^{(i+1)}y_j}{\sum_{j=1}^n\mu_j^{(i+1)}}\\ &q^{(i+1)}=\frac{\sum_{j=1}^n(1-\mu_j^{(i+1)})y_j}{\sum_{j=1}^n(1-\mu^{(i+1)})} \end{aligned}

:::

在用条件期望\mu_j代替z_j后我们的式子由隐式偶合变为显式加权，避免了隐式求解时的组合爆炸问题，如此，我们终于可以使用极大似然估计的基本思想，最大化参数，得到参数的最大似然估计了

上面问题的三个参数的最大似然估计，现在代理为对Q，也就是条件期望的最大似然估计，上面的三个参数的结果是对E步中的Q(\theta,\theta^{(i)})中的三个参数分别求导等于0得到的

所以根据y_j分别等于1和时求出对应的\mu_j，使用题中的观测数据，代入上面三个公式即可求得\pi=0.5,p=q=0.6

算法

下面正式给出EM算法
输入：
观测变量数据Y，隐变量数据Z，联合分布P(Y,Z|\theta)，条件分布P(Z|Y,\theta)
输出：
模型参数\theta
（1）选择合适的参数初值\theta^{(0)}，开始迭代
注：初值可以任选，但注意EM算法对初值是敏感的
（2）E步：记\theta^{(i)}为第i次迭代参数\theta的估计值，在第i+1次迭代的E步，计算
::: align-center
Q(\theta,\theta^{(i)})=E_Z[\log P(Y,Z|\theta)|Y,\theta^{(i)}]

:::
之后使用隐变量Z在给定观测数据Y和当前参数估计\theta^{(i)}下的条件分布P(Z|Y,\theta^{(i)})（贝叶斯公式求出）代理Q中的隐变量的期望
（3）M步：求使Q(\theta,\theta^{(i)})极大化的\theta，确定i+1次迭代的参数的估计值\theta^{(i+1)}
::: align-center
\theta^{(i+1)}=\arg\max_{\theta}Q(\theta,\theta^{(i)})

:::
（4）重复第（2）步和第（3）步，直到收敛
注：Q收敛或\theta收敛

这里的Q是EM算法的核心，称为Q函数（Q Function）,下面给出它的定义：
定义（Q函数）：完全数据的对数似然函数\log P(Y,Z|\theta)在关于给定观测数据Y和当前参数\theta^{(i)}下对未观测数据Z的条件分布P(Z|Y,\theta^{(i)})的期望称为Q函数，即
::: align-center
Q(\theta,\theta^{(i)})=E_Z[\log P(Y,Z|\theta)|Y,\theta^{(i)}]

:::

EM算法的导出

上面叙述了EM算法，它能近似实现对观测数据的极大似然估计，那么为什么呢？下面通过近似求解观测数据的对数似然函数的极大化问题来导出EM算法，由此可以清楚的看出EM算法的作用

我们面对一个含有隐变量的概率模型，目标是极大化观测数据（不完全数据）Y的对数似然函数，即极大化
::: align-center
\begin{aligned}&L(\theta)=\log P(Y|\theta)\\ &=\log(\sum_Z P(Y,Z|\theta))\\ &=\log(\sum_ZP(Y|Z,\theta)P(Z|\theta)) \end{aligned}

:::
注意到上式极大化的困难是含有未观测数据并有包含和（或积分）的对数

实际上，EM算法是通过迭代逐步近似极大化L(\theta)的，假设在i次迭代后\theta的估计值是\theta^{(i)}。我们希望新估计值\theta能使L(\theta)增加，即L(\theta)>L(\theta^{(i)})，并逐步达到极大值，为此，考虑两者的差
::: align-center
L(\theta)-L(\theta^{(i)})=\log(\sum_ZP(Y|Z,\theta)P(Z|\theta))-\log P(Y|\theta^{(i)})

:::
利用凸函数的Jessen不等式（广义形式）f(\sum_j \lambda_iy_j)\le \sum_{j}\lambda_jf(y_j)，对于凹函数\log x，不等号反向，于是得到上面差的下界
::: align-center
\begin{aligned} &L(\theta)-L(\theta^{(i)})=\log(\sum_ZP(Z|Y,\theta^{(i)})\frac{P(Y|Z,\theta)P(Z,\theta)}{P(Z|Y,\theta^{(i)})})-\log P(Y|\theta^{(i)})\\ &\ge \sum_Z P(Z|Y,\theta^{(i)})\log\frac{P(Y|Z,\theta)P(Z,\theta)}{P(Z|Y,\theta^{(i)})}-\log P(Y|\theta^{(i)})\\ &=\sum_Z P(Z|Y,\theta^{(i)})\log\frac{P(Y|Z,\theta)P(Z,\theta)}{P(Z|Y,\theta^{(i)})P(Y|\theta^{(i)})}\\ \end{aligned}

:::
令
::: align-center
B(\theta,\theta^{(i)}) =L(\theta^{(i)})+\sum_Z P(Z|Y,\theta^{(i)})\log\frac{P(Y|Z,\theta)P(Z,\theta)}{P(Z|Y,\theta^{(i)})P(Y|\theta^{(i)})}

:::
则
::: align-center
L(\theta)\ge B(\theta,\theta^{(i)})

:::
即B(\theta,\theta^{(i)})是L(\theta)的一个下界，而且由上式知
::: align-center
L(\theta^{(i)})=B(\theta^{(i)},\theta^{(i)})

:::
因此，任何可以使B(\theta,\theta^{(i)})增大的\theta，也可以使L(\theta)增大。为了使L(\theta)有尽可能大的增长，选择\theta^{(i+1)}使B(\theta,\theta^{(i)})达到极大，即
::: align-center
\begin{aligned} \theta^{(i+1)}=\arg\max_{\theta}B(\theta,\theta^{(i)}) \end{aligned}

:::
现求\theta^{(i+1)}的表达式，省去对\theta的极大化而言是常数的项，有
::: align-center
\begin{aligned} &\theta^{(i+1)}=\arg\max_{\theta}(L(\theta^{(i)})+\sum_ZP(Z|Y,\theta^{(i)}))\log\frac{P(Y|Z,\theta) P(Z|\theta)}{P(Z|Y,\theta^{(i)})P(Y|\theta^{(i)})})\\ &=\arg\max_{\theta}(\sum_Z P(Z|Y,\theta^{(i)})\log P(Y|Z,\theta) P(Z|\theta))\\ &=\arg\max_{\theta}(\sum_{Z}P(Z|Y,\theta^{(i)})\log P(Y,Z|\theta))\\ &=\arg\max_{\theta}Q(\theta,\theta^{(i)}) \end{aligned}

:::
即我们极大化\theta可以通过极大化Q(\theta,\theta^{(i)})，EM算法就是通过不断求解下界极大化逼近求解对数似然函数极大化的算法

笔者自己的理解

就笔者的理解而言，EM算法是一种估计-验证-修正的算法

例如在三硬币问题中，若我们不考虑严格的建模并计算
我们只有最后硬币正反面的数据而不知道这些硬币是B还是C，需要估计的参数是A正面朝上的概率\pi（选择硬币B的概率，同时知道正面朝上的概率也就知道朝下的概率了），B和C正面朝上的概率p/q（知道正面朝上的概率反面亦然）

EM算法考虑
（1）先设出三个参数的初值，我们不妨令它们分别为P(Z=1),P(Y=1|Z=1),P(Y=1|Z=0)，其中Z为A硬币的正/反（1/0），Y为观测到硬币正/反（1/0）
（2）之后我们用观测数据去验证Z的后验概率，例如用观测数据1验证硬币为B的后验概率为P(Z=1|Y=1)，这里我们需要用到贝叶斯公式
（3）然后我们使用P(Z=1|Y=1),P(Z=0|Y=1)和观测数据去“修正”我们我们设出的初值，具体地来说，我们考虑使用隐变量的后验概率“均值”来“逼近”真实概率
例如在修正P(Z=1)时，我们使用P(Z=1|Y)关于所有观测数据的均值，在修正P(Y=1|Z=1)时我们使用P(Z=1|Y=1)关于Y=1在P(Z=1|Y)关于所有观测数据上的均值

依笔者理解，EM算法利用的是在样本数量足够大时，后验概率均值的期望会渐趋近于真实的隐变量概率的分布，这就好比我们要把1单位的糖果分配给A，B，我们在根据他们二人对每次分到糖果的满意程度不断调整分配比例，最终我们会不断逼近理想分配比例

EM算法的收敛性

暂无计划整理

EM算法在高斯混合模型学习中的应用

高斯混合模型

高斯混合模型是指具有如下形式的概率分布模型
::: align-center
P(y|\theta)=\sum_{k=1}^K\alpha_k\phi(y|\theta_k)

:::
其中\alpha_k是系数，\alpha_k\ge 0，\sum_{k=1}^K\alpha_k=1
\phi(y|\theta_k)是高斯分布密度，\theta_k=(\mu_k,\sigma_k^2)，
::: align-center
\phi(y|\theta_k)=\frac{1}{\sqrt{2\pi}\sigma_k}\exp(-\frac{(y-\mu_k)^2}{2\sigma_k^2})

:::
称为第k个分模型

高斯混合模型参数估计的EM算法

假设观测数据y_1,y_2,\cdots,y_N由高斯混合模型生成，
::: align-center
P(y|\theta)=\sum_{k=1}^K\alpha_k\phi(y|\theta_k)

:::

其中，\theta=(\alpha_1,\alpha_2,\cdots,a_K;\theta_1,\theta_2,\cdots,\theta_K)。我们用EM算法估计高斯混合模型的参数\theta

1.明确隐变量

可以设想观测数据y_j,j=1,2,\cdots,N，是这样产生的：首先依照概率\alpha_k选择k个高斯分布分模型\phi(y|\theta_k)，然后依照第k个分模型的概率分布生成观测数据。这时观测数据y_j是已知的，反映观测数据y_j来自来自第k个分模型的数据是未知的，k=1,2,\cdots,K，以隐变量\gamma_{jk}表示
::: align-center
\gamma_{jk}=\left\{\begin{matrix} 1 & The\space j^{th}\space sample\space comes \space from\space k^{th}\space model\\ 0 & Otherwise \end{matrix}\right .

:::
\gamma_{jk}是0-1随机变量

有了观测数据y_j及未观测数据\gamma_{jk}，那么完全数据是
::: align-center
(y_j,\gamma_{j1},\gamma_{j2},\cdots,\gamma_{jK}),j=1,2,\cdots,N

:::
于是可以写出完全数据的似然函数：
::: align-center
\begin{aligned} &P(y,\gamma|\theta)=\prod_{j=1}^N P(y_j,\gamma_{j1},\gamma_{j2},\cdots,\gamma_{jK}|\theta)\\ &=\prod_{k=1}^K\prod_{j=1}^N[\alpha_k\phi(y_j|\theta_k)]^{\gamma_{jk}}\\ &=\prod_{k=1}^K\alpha^{n_k}\prod_{j=1}^N[\phi(y_j|\theta_k)]^{\gamma_{jk}}\\ &=\prod_{k=1}^K\alpha^{n_k}\prod_{i=1}^N[\frac{1}{\sqrt{2\pi}\sigma_k}\exp(-\frac{(y_j-\mu_k)^2}{2\sigma_k^2})]^{\gamma_{jk}} \end{aligned}

:::
式中，n_k=\sum_{j=1}^N\gamma_{jk},\sum_{k=1}^Kn_k=N
那么，完全数据的对数似然函数为
::: align-center
\log P(y,\gamma|\theta)=\sum_{k=1}^K\{n_k\log\alpha_k+\sum_{i=1}^N\gamma_{jk}[\log\frac{1}{\sqrt{2\pi}}-\log\sigma_k-\frac{(y_j-\mu_k)^2}{2\sigma_k^2}]\}

:::

2.EM算法的E步：确定Q函数

对对数似然函数取条件期望
::: align-center
\begin{aligned} &Q(\theta|\theta^{(t)})=E[\log P(y,\gamma|\theta)=\sum_{k=1}^K\{n_k\log\alpha_k+\sum_{i=1}^N\gamma_{jk}[\log\frac{1}{\sqrt{2\pi}}-\log\sigma_k-\frac{(y_j-\mu_k)^2}{2\sigma_k^2}]\}\\ &=\sum_{k=1}^K\{\sum_{j=1}^N(E_{\gamma_{jk}})\log\alpha_k+\sum_{j=1}^N(E_{\gamma_{jk}})[\log\frac{1}{\sqrt{2\pi}}-\log\sigma_k-\frac{(y_j-\mu_k)^2}{2\sigma_k^2}]\}]\} \end{aligned}

:::
这里需要计算条件期望E(\gamma_{jk}|y,\theta)，记为\hat{r}_{jk}
::: align-center
\begin{aligned} &\hat{r}_{jk}=P(\gamma_{jk}=1|y,\theta)\\ &=\frac{P(y_j|\gamma_{jk=1},\theta)P(\gamma_{jk}=1|\theta)}{\sum_{k=1}^K P(y_j|\gamma_{jk=1},\theta)P(\gamma_{jk}=1|\theta)}\\ &=\frac{\alpha_k\phi(y_j|\theta_k)}{\sum_{k=1}^K\alpha_k\phi(y_j|\theta_k)},j=1,2,\cdots,N;k=1,2,\cdots,K \end{aligned}

:::
\hat{\gamma}_{jk}是当前模型参数下第j个观测数据来自第k个分模型的概率，称为分模型k对观测数据y_j的响应度

将\hat{\gamma}_{jk}=E\gamma_{jk}及n_k=\sum_{j=1}^N E\gamma_{jk}代入上式得
::: align-center
Q(\theta|\theta^{(t)})=\sum_{k=1}^K\{n_k\log\alpha_k+\sum_{j=1}^N \hat{\gamma}_{jk}[\log\frac{1}{\sqrt{2\pi}}-\log\sigma_k-\frac{(y_j-\mu_k)^2}{2\sigma_k^2}]\}]\}

:::

3.确定EM算法的M步

迭代的M步是求函数Q(\theta|\theta^{(t)})对\theta的极大值，即求新一轮迭代的模型参数：
::: align-center
\theta^{(i+1)}=\arg\max_{\theta}Q(\theta|\theta^{(i)})

:::
用\hat{\mu}_k,\hat{\sigma}_{k}^2,\hat{\alpha}_k,k=1,2,\cdots,K表示\theta^{(i+1)}的各参数，求\hat{\mu}_k,\hat{\sigma}_{k}^2只需对上面的Q函数分别求偏导并令其为0，而\hat{\alpha}_k是在\sum_{k=1}^K\alpha_k=1条件下求偏导并令其为0得到的，结果如下
::: align-center
\begin{aligned} &\hat{\mu}_{k}=\frac{\sum_{j=1}^N \hat{\gamma}_{jk}y_j}{\sum_{j=1}^N \hat{\gamma}_{jk}},k=1,2,\cdots,K\\ &\sigma_k^2=\frac{\sum_{j=1}^N \hat{\gamma}_{jk}(y_j-\mu_k)^2}{\sum_{j=1}^N \hat{\gamma}_{jk}},k=1,2,\cdots,K\\ &\alpha_k=\frac{n_k}{N}=\frac{\sum_{j=1}^N \hat{r}_{jk}}{N},k=1,2,\cdots,K \end{aligned}

:::
重复以上计算，直到对数似然函数值不再有明显的变化为止。

总结

现将估计高斯混合模型参数的EM算法总结如下：
输入：观测数据y_1,y_2,\cdots,y_N，高斯混合模型
输出：高斯混合模型参数
（1）取参数的初始值开始迭代
（2）E步：依据当前模型参数，计算分模型k对观测数据y_j的响应度
::: align-center
\frac{\alpha_k\phi(y_j|\theta_k)}{\sum_{k=1}^K\alpha_k\phi(y_j|\theta_k)},j=1,2,\cdots,N;k=1,2,\cdots,K

:::
（3）M步：计算新一轮迭代的模型参数
::: align-center
\begin{aligned} &\hat{\mu}_{k}=\frac{\sum_{j=1}^N \hat{\gamma}_{jk}y_j}{\sum_{j=1}^N \hat{\gamma}_{jk}},k=1,2,\cdots,K\\ &\sigma_k^2=\frac{\sum_{j=1}^N \hat{\gamma}_{jk}(y_j-\mu_k)^2}{\sum_{j=1}^N \hat{\gamma}_{jk}},k=1,2,\cdots,K\\ &\alpha_k=\frac{n_k}{N}=\frac{\sum_{j=1}^N \hat{r}_{jk}}{N},k=1,2,\cdots,K \end{aligned}

:::
（4）重复第2步和第3步，直到收敛

EM算法的推广

暂无计划整理

简单练习

两硬币问题：设有两枚硬币A和B，正面概率\theta_a和\theta_b（未知）。每次实验随机选择一枚硬币投掷M次，重复N次实验。仅观测到每次实现的正反面次数，而不知道选择的是哪枚硬币。
目标：通过EM算法估计\theta_a,\theta_b

解：
设隐变量z_i表示第i次实验选择的硬币，z_i=1,2分别表示选到了A/B硬币
设\theta_a=0.5,\theta_b=0.5

E步：计算每次实验i选择硬币A的后验概率\omega_i：
::: align-center
\omega_i=\frac{\theta_a^{h_i}(1-\theta_a)^{t_i}}{\theta_a^{h_i}(1-\theta_a)^{t_i}+\theta_b^{h_i}(1-\theta_b)^{t_i}}

:::
h_i,t_i分别代表
M步：使用后验概率修正
::: align-center
\begin{aligned} &\theta_a^{new}=\frac{\sum_{i=1}^N \omega_ih_i}{\sum_{i=1}^N \omega_i M}\\ &\theta_b^{new}=\frac{\sum_{i=1}^N(1-\omega_i)h_i}{\sum_{i=1}^N(1-\omega_i)M} \end{aligned}

:::
具体来说：
定义Q函数
首先写出完全数据的联合概率似然函数
::: align-center
P(X,Z|\theta)=\prod_{i=1}^NP(Z_i)\cdot P(h_i,t_i|Z_i,\theta)

:::
之后取对数似然函数，并取条件期望得Q函数
::: align-center
Q(\theta|\theta^{(t)})=\sum_{i=1}^NE_{Z|X,\theta^{(t)}}[\log P(Z_i)+\log P(h_i,t_i|Z_i,\theta)]

:::
由于P(Z_i)是常数，并不影响Q函数的极大化，所以我们只考虑后一部分，令它为新的Q函数，并将E步中计算的Z_i的后验概率概率代入后展开期望为
::: align-center
\omega_i\log P(h_i,t_i|Z_i=A,\theta)+(1-\omega_i)\log P(h_i,t_i|Z_i=B,\theta)

:::
其中，P(h_i,t_i|Z_i=A)的单次实验的似然为
::: align-center
P(h_i,t_i|Z_i=A)=\theta_a^{h_i}(1-\theta_a)^{t_i}

:::
P(h_i,t_i|Z_i=B)类似

代入新的Q函数
::: align-center
\begin{aligned} &Q(\theta|\theta^{(t)})=\sum_{i=1}^N[\omega_i(h_i\log \theta_a+t_i\log(1-\theta_a))\\ &+(1-\omega_i)(h_i\log \theta_b+t_i\log(1-\theta_b))] \end{aligned}

:::
分别对\theta_a,\theta_b求偏导即可得到结果