混合高斯分布与 EM 算法

概率论中的 Jensen 不等式

对于 Jensen 不等式,通常情况下是这样的:对于 \(f^{\prime \prime}(x) \geq 0\) 也就是对于凸函数而言,这个可以用中值定理来证明,同时这个公式还可以用于证明多项式的均值不等式与调和不等式,这里就不给出证明了,在概率的条件下就是:实际上,这就是均值不等式来的啊, E 表示数学期望

\[\mathrm{E}[f(X)] \geq f(\mathrm{E} X) \]

EM 算法讲了什么

对于下面的第一行的式子, 本质上这个式子是函数隐变量的, 将这些隐变量显示的写出来就是下面的第二行的式子:

\[\begin{aligned} \ell(\theta) &=\sum_{i=1}^{N} \log p(x ; \theta) \\ &=\sum_{i=1}^{N} \log \sum_{z} p(x, z ; \theta) \end{aligned} \]

现在,我们假设对于每一个 $ i $ , \(Q_{i}\) 表示 $ z $ 的分布,那么我们有 \(\sum_{z} Q_{i}(z)=1, Q_{i}(z) \geq0\) 。然后我们使用 Jensen 不等式将上面的式子进行放缩,写成下面的这样形式,

\[\begin{aligned} \ell(\theta) =\sum_{i} \log p\left(x^{(i)} ; \theta\right) &=\sum_{i} \log \sum_{z^{(i)}} p\left(x^{(i)}, z^{(i)} ; \theta\right) \\ &=\sum_{i} \log \sum_{z^{(i)}} Q_{i}\left(z^{(i)}\right) \frac{p\left(x^{(i)}, z^{(i)} ; \theta\right)}{Q_{i}\left(z^{(i)}\right)} \\ & \geq \sum_{i} \sum_{z^{(i)}} Q_{i}\left(z^{(i)}\right) \log \frac{p\left(x^{(i)}, z^{(i)} ; \theta\right)}{Q_{i}\left(z^{(i)}\right)} \end{aligned} \]

这里的 \(Q_{i}\) 应该就是 EM 算法不断的迭代的关键。

上面的式子表明 \(\ell(\theta)\) 有一个下界,从极大似然的角度考虑,我们就是要最大化这个下界,得到 $ \theta $ 的取值,但是其中的 \(Q_{i}\) 是一个隐变量的分布,我们并不知道这个分布是什么样子的。

上面使用 Jensen 不等式关键的一步是:

\[f\left(\mathrm{E}_{z^{(i)} \sim Q_{i}}\left[\frac{p\left(x^{(i)}, z^{(i)} ; \theta\right)}{Q_{i}\left(z^{(i)}\right)}\right]\right) \geq \mathrm{E}_{z^{(i)} \sim Q_{i}}\left[f\left(\frac{p\left(x^{(i)}, z^{(i)} ; \theta\right)}{Q_{i}\left(z^{(i)}\right)}\right)\right] \]

这个式子取等式的时候,对于数学期望而言,只有常数可以满足数学期望中的 Jensen 不等式相等。这里不做具体的证明,我们可以考虑从均值不等式来理解这个问题,假设这个常数是 $ c$ :

\[\frac{p\left(x^{(i)}, z^{(i)} ; \theta\right)}{Q_{i}\left(z^{(i)}\right)}=c \]

也就是说:

\[Q_{i}\left(z^{(i)}\right) \propto p\left(x^{(i)}, z^{(i)} ; \theta\right) \]

我们知道:\(\sum_{z} Q_{i}\left(z^{(i)}\right)=1\) ,那么我可以考虑这样的情况:

\[1 = \sum{Q_i\left( z^{\left( i \right)} \right)} \propto \sum{p\left( x^{\left( i \right)},z^{\left( i \right)};\theta \right) } = \sum_{z} p\left(x^{(i)}, z^{\left(i \right)} ;\ \theta\right) \]

这样的话,就有下面的公式:

\[\begin{aligned} Q_{i}\left(z^{(i)}\right) &=\frac{p\left(x^{(i)}, z^{(i)} ; \theta\right)}{\sum_{z} p\left(x^{(i)}, z^{\left(i \right)} ; \theta\right)} \\ &=\frac{p\left(x^{(i)}, z^{(i)} ; \theta\right)}{p\left(x^{(i)} ; \theta\right)} \\ &=p\left(z^{(i)} | x^{(i)} ; \theta\right) \end{aligned} \]

这一步我们称为 E 步骤,可以得到等式成立的条件, \(Q_{i}\left(z^{(i)}\right) :=p\left(z^{(i)} | x^{(i)} ; \theta\right)\) ,下界成立的条件,我们将这个条件带入 $\ell(\theta) $ 得到的结果就是下界,也就是说,这个是最大似然函数的条件之一, 我们需要最大似然的就是:那么我们可以用下面的步骤来计算 $ \theta$ ,

\[\theta :=\arg \max _{\theta} \sum_{i} \sum_{z^{(i)}} Q_{i}\left(z^{(i)}\right) \log \frac{p\left(x^{(i)}, z^{(i)} ; \theta\right)}{Q_{i}\left(z^{(i)}\right)} \]

那么这个 $ \theta$ 是不是最优的呢?接下来我们来证明 EM 算法中一个很关键的问题,那就是上述的 E 步骤与 M 步骤不断放缩与收敛,最终得到最优的 $ \theta$ 的值。

所以 EM 算法的步骤可以表示成下面这样:

  1. 函数声明 \(\operatorname{EM}\left(p_{Y, C}(y, c ; \theta), p_{C | Y}(c | y ; \theta), \theta^{(0)}\right)\)
  2. for iteration \(t \in 1,2, \ldots\) do
  3. \(Q_{i}^{(t)} \leftarrow P\left(z_i | x_i ; \theta^{(t-1)}\right) \quad(\text { E-step })\)
  4. \(\theta^{(t)} \leftarrow \operatorname{argmax}_{\theta} \mathbb{E}_{Q_{i}^{(t)}}\left[P(y, C ; \theta)\right] \quad\left(\mathrm{M}{-\mathrm{Step}}\right)\)
  5. if \(\theta^{(t)} \approx \theta^{(t-1)}\) then
  6. return $\theta^{(t)} $

EM 算法的收敛问题证明

这个收敛的思想是这样的:我们的 $ \theta $ 是直接通过极大似然函数算出的来的。那么 EM 算法迭代的步骤就是,我们不断地最大化极大似然估计函数,也就是说 我们每次都最大化了 \(\ell\left(\theta^{(t)}\right)\) 的下界,只有保证 \(\ell\left(\theta^{(t)}\right) \leq \ell\left(\theta^{(t+1)}\right)\) , 这样就会不断地逼近最优解, 说明算法是正确的。在 EM 迭代的过程中,我们不断改变是 \(Q_{i}\) 也就是分布函数的选择(我们本质是更新 \(\theta\) , 改变 \(\theta\) 等价于改变 \(Q_i\) ),

为了更好的说明,假设上一步我们已经得到了一个最优解 \(\ell\left(\theta^{(t)}\right)\) 以及 $ \theta^{(t)} $ 那么在这一步,我们满足下面的情况:

\[Q_{i}^{(t)}\left(z^{(i)}\right) :=p\left(z^{(i)} | x^{(i)} ; \theta^{(t)}\right) \ \ \ Jensen不等式条件\\ 用于计算 \ell\left(\theta^{(t+1)}\right) \]

\[\ell\left(\theta^{(t)}\right)=\sum_{i} \sum_{z^{(i)}} Q_{i}^{(t)}\left(z^{(i)}\right) \log \frac{p\left(x^{(i)}, z^{(i)} ; \theta^{(t)}\right)}{Q_{i}^{(t)}\left(z^{(i)}\right)} \]

而更新的参数 \(\theta^{(t+1)}\) 就是来自最大化上面右边的式子, 那么我们可以推出下面的不等式:

\[\begin{aligned} \ell\left(\theta^{(t+1)}\right) & \geq \sum_{i} \sum_{z^{(i)}} Q_{i}^{(t)}\left(z^{(i)}\right) \log \frac{p\left(x^{(i)}, z^{(i)} ; \theta^{(t+1)}\right)}{Q_{i}^{(t)}\left(z^{(i)}\right)} \\ & \geq \sum_{i} \sum_{z^{(i)}} Q_{i}^{(t)}\left(z^{(i)}\right) \log \frac{p\left(x^{(i)}, z^{(i)} ; \theta^{(t)}\right)}{Q_{i}^{(t)}\left(z^{(i)}\right)} \\ &=\ell\left(\theta^{(t)}\right) \end{aligned} \]

对于上面的两个不等式做出以下解释:

  1. 第一个不等式来自 Jensen 不等式,原式是这样的:\(\ell(\theta) \geq \sum_{i} \sum_{z^{(i)}} Q_{i}\left(z^{(i)}\right) \log \frac{p\left(x^{(i)}, z^{(i)} ; \theta\right)}{Q_{i}\left(z^{(i)}\right)}\) 对任意的$ Q_i $ 与\(\theta\) 均成立,这里为什么对任意的 \(\theta\) 都成立, 是因为 \(Q\) 本质也是由 \(\theta\) 决定的.
  2. 对于第二个不等式, 就是我们前面 M 步骤的结果, 最大化下面的极大似然函数的时候得到$ \theta^{(t+1)}$ 所以第二个式子显然:

\[\sum_{i} \sum_{z^{(i)}} Q_{i}^{(t)}\left(z^{(i)}\right) \log \frac{p\left(x^{(i)}, z^{(i)} ; \theta^{(t)}\right)}{Q_{i}^{(t)}\left(z^{(i)}\right)} \]

综上证明了 EM 算法的收敛性问题

混合高斯分布

在一般的情况下,对于所得到的样本集,\(X=\left\{x_{1}, \dots, x_{N}\right\}\),我们的目标是最大化似然函数,通过最大化似然函数来获取参数的值。这是似然函数往往取对数表示就是:

\[\begin{aligned} L(\theta | X) &=\log \left(\prod_{i=1}^{N} p\left(x_{i} | \theta\right)\right) \\ &=\sum_{i=1}^{N} \log p\left(x_{i} | \theta\right) \end{aligned} \]

这个过程的参数估计可以描述成:

\[\hat{\theta}=\arg \max _{\theta} L(\theta | X) \]

这个结果是可以直接计算出来的,那么引入混合高斯分布会是什么样呢?

混合高斯分布:

简单的说就是,非混合的情况下,我们的数据集满足高斯分布,这样用极大似然就直接算出高斯分布中的参数就可以了。那么混合的情况下就是,数据集是由多个满足高斯分布的数据集线性组合起来,这个时候我们理解为:有 \(k\) 个满足不同参数的高斯分布的原数据集,并且 \(\sum_{j=1}^{k} \pi_{j}=1\).

此时:

\[p(\boldsymbol{x})=\sum_{k=1}^{K} \pi_{k} \mathcal{N}\left(\boldsymbol{x} | \boldsymbol{\mu}_{k}, \boldsymbol{\Sigma}_{k}\right) \]

举个例子说明:

我们假设男生的身高符合正态分布, \(\mathcal{N_1}\left(\boldsymbol{\mu}_{1}, \boldsymbol{\Sigma}_{1}\right)\) 女生的身高符合正态分布 \(\mathcal{N_2}\left(\boldsymbol{\mu}_{2}, \boldsymbol{\Sigma}_{2}\right)\) , 现在我们得到一个关于身高的样本数据集, 但是不知道男生与女生的比例, 要拟合这个身高数据集符合哪种分布, 这里就可以使用混合高斯分布,

\[p(\boldsymbol{x})= \pi_{1} \mathcal{N}\left(\boldsymbol{x} | \boldsymbol{\mu}_{1}, \boldsymbol{\Sigma}_{1}\right) + \pi_{2} \mathcal{N}\left(\boldsymbol{x} | \boldsymbol{\mu}_{2}, \boldsymbol{\Sigma}_{2}\right) \]

其中 \(\pi_1\)\(\pi_2\) 就表示男生与女生的比例. 也表示 \(p(x)\) 对不同高斯分布的可能性,

在上面的问题中, \(\pi_1\)\(\pi_2\) 是未知参数 , 那么这个例子与一个问题很类似, 那就是贝叶斯分布的问题, 如果我们在已知后验概率的情况下, 就可以很简单的计算先验概率了, 也就是说, 如果我们已知男生满足的高斯分布的参数, \(\mu_1\)\({\Sigma}_{1}\) 和 女生高斯分布的参数 \(\mu_2\)\({\Sigma}_{2}\) , 这时就可以反过来使用贝叶斯计算出 \(\pi_1\)\(\pi_2\) 的值.

混合高斯分布与贝叶斯

对于一个混合高斯分布来说, \(\pi_1, \pi_2, \dots, \pi_n\) 对这个分布来说也是未知参数, 这里我们从贝叶斯的角度来改写这个未知参数的形式, 这样做的目的是为了更好的使用下面的 EM 算法以及 EM 算法的计算方式, 对于上面的 \(N\) 维混合高斯分布, 我们取一个向量\(N\) 维向量 \({z}\) 表示取得的, 其中 $ z_k(1 \leq k \leq K) $, 只能取 0 与 1 两个值, 其中 \(p(z_i == 1)\) 表示第 i 中高斯分布被选中的概率, 而 \(p(z_i = 0)\) 表示未被选中的概率, 那么用向量 \(z\) 表示混合高斯分布就有两种方式, 一种是类似于

\[[1,1,1,0,0,\dots, 1, 0, 1] \]

的形式, 这样直接使用向量的乘法的到概率,

而另一种是

\[[0,1,0, 0\dots, \dots, 0] \]

的形式, 这里只有一个 \(z_i\) 为 1, 因为我们要使用贝叶斯分布, 反过来使用全概率公式, 所以我们使用第二种形式, 在这种形式下, \(z_i\) 之间是相互独立的, 所以 \(z\) 向量的概率可以写成:

\[p(\boldsymbol{z})=p\left(z_{1}\right) p\left(z_{2}\right) \ldots p\left(z_{K}\right)=\prod_{k=1}^{K} \pi_{k}^{z_{k}} \]

然后对应于 \(z_i\) 的那一个分布是服从高斯分布的, 可以将对应的 \(z_i\) 分布的概率表示如下:

\[p\left(\boldsymbol{x} | z_{i}=1\right)=\mathcal{N}\left(\boldsymbol{x} | \boldsymbol{\mu}_{i}, \boldsymbol{\Sigma}_{i}\right) \]

进而上式有可以写成如下形式:

\[p(\boldsymbol{x} | \boldsymbol{z})=\prod_{i=1}^{N} \mathcal{N}\left(\boldsymbol{x} | \boldsymbol{\mu}_{i}, \boldsymbol{\Sigma}_{i}\right)^{z_{i}} \]

我们知道, 在向量 \(z\) 中, 假设 \(z_j = 1\) 那么对于所有的 \(z_i (i\neq j)\) 都有 \(z_i = 0\), 所以上式成立很显然, 我们得到了 \(p(z)\)\(p(x|z)\) 就可以得到 \(p(x)\) 的公式了:

\[\begin{aligned} p(\boldsymbol{x}) &=\sum_{\boldsymbol{z}} p(\boldsymbol{z}) p(\boldsymbol{x} | \boldsymbol{z}) \\ &=\sum_{\boldsymbol{z}}\left(\prod_{k=1}^{K} \pi_{k}^{z_{k}} \mathcal{N}\left(\boldsymbol{x} | \boldsymbol{\mu}_{k}, \boldsymbol{\Sigma}_{k}\right)^{z_{k}}\right) \\ &=\sum_{k=1}^{K} \pi_{k} \mathcal{N}\left(\boldsymbol{x} | \boldsymbol{\mu}_{k}, \boldsymbol{\Sigma}_{k}\right) \end{aligned} \tag{1} \]

公式 (1) 本质上使用的是全概率公式, 因为向量 \(z\) 的表示方式我们取的是第二种, 上式中 \(p(z)p(x|z)\) 的结果是 \(p(x,z)\) 所以要用全概率公式相加, 从上面的式子我们惊奇的发现, 我们 \(p(x)\) 改写成了后验概率的形式, 因为我们划分了所有的 \(z_i\) , 在不同的向量 \(z\) 的基础上得到的概率 \(p(x)\), 这时反过来就可以求先验概率:

\[\begin{aligned} \gamma\left(z_{i}\right) &=p\left(z_{i}=1 | \boldsymbol{x}\right) \\ &=\frac{p\left(z_{i}=1\right) p\left(\boldsymbol{x} | z_{i}=1\right)}{p\left(\boldsymbol{x}, z_{i}=1\right)} \\ &=\frac{p\left(z_{i}=1\right) p\left(\boldsymbol{x} | z_{i}=1\right)}{\sum_{j=1}^{N} p\left(z_{j}=1\right) p\left(\boldsymbol{x} | z_{j}=1\right)}(\text { 全概率公式 }) \\ &=\frac{\pi_{i} \mathcal{N}\left(\boldsymbol{x} | \boldsymbol{\mu}_{i}, \boldsymbol{\Sigma}_{\boldsymbol{i}}\right)}{\sum_{j=1}^{N} \pi_{j} \mathcal{N}\left(\boldsymbol{x} | \boldsymbol{\mu}_{j}, \boldsymbol{\Sigma}_{j}\right)}(\text {结合}(1) . \end{aligned} \tag{2} \]

EM 算法用于混合高斯分布

前面的 EM 算法是一般情况下的 EM 算法求解以及证明的过程, 那么在实际的计算的过程中, EM 算法求解 GMM 问题是什么样的呢? 下面的公式 (3) 是基本的公式,

\[p(\boldsymbol{x|\pi, \mu, \Sigma})=\sum_{k=1}^{K} \pi_{i} \mathcal{N}\left(\boldsymbol{x} | \boldsymbol{\mu}_{i}, \boldsymbol{\Sigma}_{i}\right) \tag{3} \]

这个公式与 EM 算法最终的最大化的公式是一样的(因为 Jensen 不等式的条件), 所以对这个公式求极大似然,

\[\prod_{i}^{N} p(\boldsymbol{x|\pi, \mu, \Sigma}) = \prod_{n}^{N} \sum_{i}^{N} \pi_{i} \mathcal{N}\left(\boldsymbol{x} | \boldsymbol{\mu}_{i}, \boldsymbol{\Sigma}_{i}\right) \tag{4} \]

然后对这个公式求对数, 以及求偏导, 可以很容易得到下面的公式:

\[0=-\sum_{n=1}^{N} \frac{\pi_{i} \mathcal{N}\left(\boldsymbol{x}_{n} | \boldsymbol{\mu}_{i}, \boldsymbol{\Sigma}_{i}\right)}{\sum_{j} \pi_{j} \mathcal{N}\left(\boldsymbol{x}_{n} | \boldsymbol{\mu}_{j}, \boldsymbol{\Sigma}_{j}\right)} \boldsymbol{\Sigma}_{i}^{-1}\left(\boldsymbol{x}_{n}-\boldsymbol{\mu}_{i}\right) \tag{5} \]

对于公式(5), 前面部分是通常的对数求导, 将 \(\sum^{K} \pi_{i} \mathcal{N}\left(\boldsymbol{x} | \boldsymbol{\mu}_{i}, \boldsymbol{\Sigma}_{i}\right)\) 放在分母, 然后后面的 \(\boldsymbol{\Sigma}_{i}^{-1}\left(\boldsymbol{x}_{n}-\boldsymbol{\mu}_{i}\right)\) 就是一维的正态分布的求导公式, 这样就可以求出 \(\mu_i\) 的值了,

\[\boldsymbol{\mu}_{i}=\frac{1}{N_{i}} \sum_{n=1}^{N} \gamma\left(z_{n i}\right) \boldsymbol{x}_{n} \tag{6.1} \]

其中

\[N_i = \sum_{n=1}^N \gamma\left(z_{n i}\right) \tag{6.2} \]

对上面公式的解释:

  1. \(\gamma\left(z_{n i}\right)\) 表示对于样本 \(x_n\) 属于第 \(i\) 类的概率, 在男生女生问题中就是对于某一个样本属于男生或者女生的概率
  2. \(N_i\) 就是对于所有样本, 属于第 \(i\) 类的概率的个数总和
  3. \(\mu_i\) 就将这个值进行加权平均, 在公式 6.1 中, \(x_n\) 表示权值, \(\gamma\left(z_{n i}\right)\) 表示概率

同理, 使用极大似然的方法可以求得:

\[\boldsymbol{\Sigma}_{i}=\frac{1}{N_{i}} \sum_{n=1}^{N} \gamma\left(z_{n i}\right)\left(x_{n}-\boldsymbol{\mu}_{i}\right)\left(x_{n}-\boldsymbol{\mu}_{i}\right)^{T} \]

最后是 \(\pi_1, \pi_2, \dots, \pi_n\) , 这个的极大值是条件极大值, 条件限制是 \(\sum_{1}^{N} \pi_i = 1\) , 也就是在求偏导的时候需要加入拉格朗日算子, 使用拉格朗日乘子法:

\[\ln p(\boldsymbol{x} | \boldsymbol{\pi}, \boldsymbol{\mu}, \boldsymbol{\Sigma})+\lambda\left(\sum_{i=1}^{N} \pi_{i}-1\right) \tag{7.1} \]

求上式关于 \(\pi_i\) 的极大似然函数, 得到:

\[0=\sum_{n=1}^{N} \frac{\mathcal{N}\left(\boldsymbol{x}_{n} | \boldsymbol{\mu}_{i}, \boldsymbol{\Sigma}_{i}\right)}{\sum_{j} \pi_{j} \mathcal{N}\left(\boldsymbol{x}_{n} | \boldsymbol{\mu}_{j}, \boldsymbol{\Sigma}_{j}\right)}+\lambda \tag{7.2} \]

公式 7.2 求导的过程就是, 第一部分与前面对 \(\mu\) 求导类似, 这里分子不变很显然, 对于结果两边同乘以 \(\pi_i\) 得到:

\[0=\sum_{n=1}^{N} \frac{\pi_{i} \mathcal{N}\left(\boldsymbol{x}_{n} | \boldsymbol{\mu}_{i}, \boldsymbol{\Sigma}_{i}\right)}{\sum_{j} \pi_{j} \mathcal{N}\left(\boldsymbol{x}_{n} | \boldsymbol{\mu}_{j}, \boldsymbol{\Sigma}_{j}\right)}+\lambda \pi_{i} \]

结合公式 (2) 可以得到:

\[0 = N_i + \lambda \pi_i \\ 0 = \sum_{i = 1}^{N} N_i + \lambda \sum_{i=1}^{N} \pi_i \\ 0 = N+k \tag{7.3} \]

最终得到就是:

\[\pi_i = \frac{N_i}{N} \]

EM 算法求解的迭代过程

EM 算法的不断迭代收敛过程就是不断地计算 E 步骤与 M 步骤, 直到期望的似然函数收敛,

E Step

根据当前的 \(\pi_i\) , \(\mu_i\) , \(\Sigma_i\) 计算:

\[y(z_nk) = \frac{\pi_{i} \mathcal{N}\left(\boldsymbol{x}_{n} | \boldsymbol{\mu}_{i}, \boldsymbol{\Sigma}_{i}\right)}{\sum_{j}^{N} \pi_{j} \mathcal{N}\left(\boldsymbol{x}_{n} | \boldsymbol{\mu}_{j}, \boldsymbol{\Sigma}_{j}\right)} \]

同时需要记录的值为:

\[p(\boldsymbol{x})=\sum_{i=1}^{N} \pi_{i} \mathcal{N}\left(\boldsymbol{x} | \boldsymbol{\mu}_{i}, \boldsymbol{\Sigma}_{i}\right) \]

M Step

使用上面的公式更新参数:

\[\boldsymbol{\mu}_{i+1}=\frac{1}{N_{i}} \sum_{n=1}^{N} \gamma\left(z_{n i}\right) \boldsymbol{x}_{n} \\ \boldsymbol{\Sigma}_{i+1}=\frac{1}{N_{i}} \sum_{n=1}^{N} \gamma\left(z_{n i}\right)\left(x_{n}-\boldsymbol{\mu}_{i}\right)\left(x_{n}-\boldsymbol{\mu}_{i}\right)^{T} \\\\ \pi_{i+1} = \frac{N_i}{N} \tag{8} \]

其中:

\[N_i = \sum_{n=1}^N \gamma\left(z_{n i}\right) \tag{6.2} \]

最后再计算,

\[\ln p(\boldsymbol{x} | \boldsymbol{\pi}, \boldsymbol{\mu}, \boldsymbol{\Sigma})=\sum_{n=1}^{N} \ln \left\{\sum_{i=1}^{N} \pi_{i} \mathcal{N}\left(\boldsymbol{x}_{i} | \boldsymbol{\mu}_{i}, \boldsymbol{\Sigma}_{i}\right)\right\} \]

判断是否收敛, 不收敛, 返回至 E Step, 然后在 M Step 继续更新, 直到收敛.

posted @ 2019-06-17 21:06  虾野百鹤  阅读(1324)  评论(0编辑  收藏  举报