问题起源
在统计学中，对于很多的样本，可以用高斯分布去描述其分布，高斯分布非常通用，但是如果现在有一组可视化的二维数据，其分布为：
很明显，当尝试使用二元高斯分布描述是不合适的，单个高斯分布无法合理描述上图的样本分布。
为了解决这个困难，引入高斯混合模型：
- 高斯指的是基础模型还是高斯分布；
- 混合指的是利用多个高斯分布加权叠加，也就是把多个不同的高斯分布的概率密度函数加权，形成一个新的概率密度函数：
p ( x ) = ∑ k = 1 K α k N ( μ k , Σ k ) , ∑ k = 1 K α k = 1 p(x)=\sum_{k=1}^{K}\alpha_{k}N(\mu_{k},\Sigma_{k}),\sum_{k=1}^{K}\alpha_{k}=1 p(x)=k=1∑KαkN(μk,Σk),k=1∑Kαk=1
假设混合模型中有
K
K
K个高斯分布，每个高斯分布的概率密度函数
N
(
μ
k
,
Σ
k
)
N(\mu_{k},\Sigma_{k})
N(μk,Σk)的权重为
α
k
\alpha_{k}
αk，针对上面的例子，可以假设混合模型有3个高斯分布，若使用不同的颜色区分不同的高斯分布，可以得到下图：
各个分布的权重为
α
1
,
α
2
,
α
3
\alpha_{1},\alpha_{2},\alpha_{3}
α1,α2,α3，将其加权叠加得到所有样本的最终概率密度。即：对于一个样本，利用高斯混合模型中的每一个高斯分布，都计算出一个概率密度值，再根据权重加权得到最终的实际概率密度值。以上是基于几何角度的分析。
从混合模型角度分析
下面从混合模型的角度，即生成模型的角度进行分析，揭示模型内部的机理：
首先，高斯混合模型中也是同时含有观测变量与隐变量的。观测变量 x x x容易理解，即样本里各个特征的观测值，在二维的可视化中，就是各个点的坐标。
而隐变量 z z z是什么？我们先用软分类的思想做一个类比，对于任意一个样本，他并不是"非此即彼"或"生硬"地属于某一个高斯分布，而是同时属于多个高斯分布，而该样本属于各个高斯分布的概率有所不同。
假设有 K K K个高斯分布，即 C 1 , C 2 , . . . , C K C_{1},C_{2},...,C_{K} C1,C2,...,CK，对于一个样本，它属于这 K K K个高斯分布的概率为 p 1 , p 2 , . . . , p K p_{1},p_{2},...,p_{K} p1,p2,...,pK，每个样本对应某个高斯分布，这就是隐变量 z z z，这个隐变量是一个离散的随机变量，其分布如下：
| z z z的取值类型 | C 1 C_{1} C1 | C 2 C_{2} C2 | … | C K C_{K} CK |
|---|---|---|---|---|
| p ( z ) p(z) p(z) | p 1 p_{1} p1 | p 2 p_{2} p2 | … | p K p_{K} pK |
现在，从混合模型的角度梳理一下生成样本的过程：
第一步：生成一个隐变量
z
z
z，换言之，依照
z
z
z的随机变量分布，依概率决定当前这个样本属于哪一个高斯分布，想象一个场景，有一个不均匀的骰子，他有
K
K
K个面，每个面向上的概率依次为
p
1
,
p
2
,
.
.
.
,
p
K
p_{1},p_{2},...,p_{K}
p1,p2,...,pK，扔一次骰子，哪个面向上，就决定样本属于哪个高斯分布，记为
C
i
C_{i}
Ci；
第二步：生成观测变量
x
x
x，此时依据决定出的高斯分布
C
i
C_{i}
Ci，生成一个样本，得到观测值；混合模型样本生成的过程看起来与前面提到概率密度加权好像没什么关系，先从混合模型角度看概率密度函数：
p
(
x
)
=
∑
z
p
(
x
,
z
)
=
∑
k
=
1
K
p
(
x
,
z
=
C
k
)
=
∑
k
=
1
K
p
(
z
=
C
k
)
p
(
x
∣
z
=
C
k
)
p(x)=\sum_{z}p(x,z)=\sum_{k=1}^{K}p(x,z=C_{k})=\sum_{k=1}^{K}p(z=C_{k})p(x|z=C_{k})
p(x)=z∑p(x,z)=k=1∑Kp(x,z=Ck)=k=1∑Kp(z=Ck)p(x∣z=Ck)
=
∑
k
=
1
K
p
k
N
(
x
∣
μ
k
,
Σ
k
)
=\sum_{k=1}^{K}p_{k}N(x|\mu_{k},\Sigma_{k})
=k=1∑KpkN(x∣μk,Σk)
此时发现，不论是从加权平均还是混合模型生成的角度理解高斯混合模型，得到的结果都是等效的。
高斯混合模型的参数估计尝试
现在已梳理了高斯混合模型的由来，所以下一阶段是参数的估计，高斯混合模型的参数包含三个部分：样本属于每个高斯分布的概率，每个高斯分布的均值向量，每个高斯分布的协方差矩阵：
θ
=
(
p
1
,
.
.
.
,
p
K
,
μ
1
,
.
.
.
,
μ
K
,
Σ
1
,
.
.
.
,
Σ
K
)
\theta=(p_{1},...,p_{K},\mu_{1},...,\mu_{K},\Sigma_{1},...,\Sigma_{K})
θ=(p1,...,pK,μ1,...,μK,Σ1,...,ΣK)
样本
X
X
X称为观测数据，
Z
Z
Z称为隐变量，
(
X
,
Z
)
(X,Z)
(X,Z)称为完全数据。
针对高斯混合模型，如果直接用极大似然法估计参数：
θ
m
l
e
=
a
r
g
m
a
x
θ
l
o
g
(
p
(
X
∣
θ
)
)
=
a
r
g
m
a
x
θ
∑
i
=
1
N
l
o
g
(
p
(
x
i
∣
θ
)
)
\theta_{mle}=argmax_{\theta}log(p(X|\theta))=argmax_{\theta}\sum_{i=1}^{N}log(p(x_{i}|\theta))
θmle=argmaxθlog(p(X∣θ))=argmaxθi=1∑Nlog(p(xi∣θ))
=
a
r
g
m
a
x
θ
∑
i
=
1
N
l
o
g
(
∑
k
=
1
K
p
k
N
(
x
i
∣
μ
k
,
Σ
k
)
)
=argmax_{\theta}\sum_{i=1}^{N}log(\sum_{k=1}^{K}p_{k}N(x_{i}|\mu_{k},\Sigma_{k}))
=argmaxθi=1∑Nlog(k=1∑KpkN(xi∣μk,Σk))
我们需要利用：
∑
i
=
1
N
l
o
g
(
∑
k
=
1
K
p
k
N
(
x
i
∣
μ
k
,
Σ
k
)
)
\sum_{i=1}^{N}log(\sum_{k=1}^{K}p_{k}N(x_{i}|\mu_{k},\Sigma_{k}))
i=1∑Nlog(k=1∑KpkN(xi∣μk,Σk))
对参数
p
k
,
μ
k
,
Σ
k
p_{k},\mu_{k},\Sigma_{k}
pk,μk,Σk分别计算偏导，但在式子中，
l
o
g
log
log中有加号，求偏导为零的根非常复杂，可以说是不能通过直接求导的方式估计参数。因此，难以使用极大似然法解偏导数的方式估计参数。解析解求不出，于是用数值解的方式，通过EM算法不断迭代。
补充：
- 高斯混合模型的参数估计可以用EM算法，也可以用梯度下降法，极大似然处理不了的原因在于求不出似然函数关于参数的偏导为0时的根，但其偏导是可以计算出的。
- 偏导数只能表示多元函数沿某个坐标轴方向的导数，梯度可以定义为一个函数的全部偏导数构成的向量，梯度向量的方向即为函数值增长最快的方向，负梯度只是给出损失函数下降最快的方向，具体下降的步长取决于学习率。
EM算法对高斯混合模型进行参数估计
下面介绍使用EM算法对高斯混合模型进行参数估计，首先定义标记：
- 样本的观测数据集： X = { x 1 , x 2 , . . . , x N } X=\left\{x_{1},x_{2},...,x_{N}\right\} X={x1,x2,...,xN}
- 隐变量数据： Z = { z 1 , z 2 , . . . , z N } Z=\left\{z_{1},z_{2},...,z_{N}\right\} Z={z1,z2,...,zN}， x i x_{i} xi与 z i z_{i} zi组成一个完整样本数据， z i z_{i} zi可取的值有 K K K类： C = { C 1 , C 2 , . . . , C K } C=\left\{C_{1},C_{2},...,C_{K}\right\} C={C1,C2,...,CK}
- 完整数据： ( X , Z ) (X,Z) (X,Z)
- 参数为： θ = ( p 1 , . . . , p K , μ 1 , . . . , μ K , Σ 1 , . . . , Σ K ) \theta=(p_{1},...,p_{K},\mu_{1},...,\mu_{K},\Sigma_{1},...,\Sigma_{K}) θ=(p1,...,pK,μ1,...,μK,Σ1,...,ΣK)
根据EM算法公式，有：
θ
t
+
1
=
a
r
g
m
a
x
θ
∫
Z
l
o
g
(
p
(
X
,
Z
∣
θ
)
)
p
(
Z
∣
X
,
θ
t
)
d
Z
\theta^{t+1}=argmax_{\theta}\int_{Z}log(p(X,Z|\theta))p(Z|X,\theta^{t})dZ
θt+1=argmaxθ∫Zlog(p(X,Z∣θ))p(Z∣X,θt)dZ
令：
Q
(
θ
,
θ
t
)
=
∫
Z
l
o
g
(
p
(
X
,
Z
∣
θ
)
)
p
(
Z
∣
X
,
θ
t
)
d
Z
Q(\theta,\theta^{t})=\int_{Z}log(p(X,Z|\theta))p(Z|X,\theta^{t})dZ
Q(θ,θt)=∫Zlog(p(X,Z∣θ))p(Z∣X,θt)dZ
由于
Z
Z
Z是离散型变量，因此转化为：
Q
(
θ
,
θ
t
)
=
∑
Z
∈
C
[
l
o
g
(
∏
i
=
1
N
p
(
x
i
,
z
i
∣
θ
)
)
⋅
∏
i
=
1
N
p
(
z
i
∣
x
i
,
θ
t
)
]
=
∑
z
i
=
C
1
,
C
2
,
.
.
.
,
C
K
[
(
∑
i
=
1
N
l
o
g
(
p
(
x
i
,
z
i
∣
θ
)
)
)
(
∏
i
=
1
N
p
(
z
i
∣
x
i
,
θ
t
)
)
]
Q(\theta,\theta^{t})=\sum_{Z\in C}[log(\prod_{i=1}^{N}p(x_{i},z_{i}|\theta))\cdot \prod_{i=1}^{N}p(z_{i}|x_{i},\theta^{t})]=\sum_{z_{i}=C_{1},C_{2},...,C_{K}}[(\sum_{i=1}^{N}log(p(x_{i},z_{i}|\theta)))(\prod_{i=1}^{N}p(z_{i}|x_{i},\theta^{t}))]
Q(θ,θt)=Z∈C∑[log(i=1∏Np(xi,zi∣θ))⋅i=1∏Np(zi∣xi,θt)]=zi=C1,C2,...,CK∑[(i=1∑Nlog(p(xi,zi∣θ)))(i=1∏Np(zi∣xi,θt))]
关于上式的说明，隐变量本身就是观测不到的，所以在计算
(
∑
i
=
1
N
l
o
g
(
p
(
x
i
,
z
i
∣
θ
)
)
)
(
∏
i
=
1
N
p
(
z
i
∣
x
i
,
θ
t
)
)
(\sum_{i=1}^{N}log(p(x_{i},z_{i}|\theta)))(\prod_{i=1}^{N}p(z_{i}|x_{i},\theta^{t}))
(∑i=1Nlog(p(xi,zi∣θ)))(∏i=1Np(zi∣xi,θt))这一项时，每次都假设
z
i
=
C
k
z_{i}=C_{k}
zi=Ck，直至
z
i
z_{i}
zi遍历完整个
C
C
C集合，再将
K
K
K个项相加。
将
Q
(
θ
,
θ
t
)
Q(\theta,\theta^{t})
Q(θ,θt)化简得到：
Q
(
θ
,
θ
t
)
=
∑
i
=
1
N
∑
z
i
=
C
1
,
.
.
.
,
C
K
[
l
o
g
(
p
(
x
i
,
z
i
∣
θ
)
)
p
(
z
i
∣
x
i
,
θ
t
)
]
Q(\theta,\theta^{t})=\sum_{i=1}^{N}\sum_{z_{i}=C_{1},...,C_{K}}[log(p(x_{i},z_{i}|\theta))p(z_{i}|x_{i},\theta^{t})]
Q(θ,θt)=i=1∑Nzi=C1,...,CK∑[log(p(xi,zi∣θ))p(zi∣xi,θt)]
此时，对数内的加号已经被转换至对数外，从而方便使用极大似然法。下面解释
p
(
x
i
,
z
i
∣
θ
)
p(x_{i},z_{i}|\theta)
p(xi,zi∣θ)和
p
(
z
i
∣
x
i
,
θ
t
)
p(z_{i}|x_{i},\theta^{t})
p(zi∣xi,θt)具体是什么：
- 首先已知：
p ( x ∣ θ ) = ∑ k = 1 K p k N ( x ∣ μ k , Σ k ) p(x|\theta)=\sum_{k=1}^{K}p_{k}N(x|\mu_{k},\Sigma_{k}) p(x∣θ)=k=1∑KpkN(x∣μk,Σk) - 对于
p
(
x
,
z
∣
θ
)
p(x,z|\theta)
p(x,z∣θ)：
p ( x , z ∣ θ ) = p ( z ∣ θ ) p ( x ∣ z , θ ) = p z N ( x ∣ μ z , Σ z ) p(x,z|\theta)=p(z|\theta)p(x|z,\theta)=p_{z}N(x|\mu_{z},\Sigma_{z}) p(x,z∣θ)=p(z∣θ)p(x∣z,θ)=pzN(x∣μz,Σz)
其中， z z z表示取第几个高斯分布； - 对于
p
(
z
∣
x
,
θ
)
p(z|x,\theta)
p(z∣x,θ)：
p ( z ∣ x , θ ) = p ( x , z ∣ θ ) p ( x ∣ θ ) = p z N ( x ∣ μ z , Σ z ) ∑ k = 1 K p k N ( x ∣ μ k , Σ k ) p(z|x,\theta)=\frac{p(x,z|\theta)}{p(x|\theta)}=\frac{p_{z}N(x|\mu_{z},\Sigma_{z})}{\sum_{k=1}^{K}p_{k}N(x|\mu_{k},\Sigma_{k})} p(z∣x,θ)=p(x∣θ)p(x,z∣θ)=∑k=1KpkN(x∣μk,Σk)pzN(x∣μz,Σz)
代入到
Q
(
θ
,
θ
t
)
Q(\theta,\theta^{t})
Q(θ,θt)为：
Q
(
θ
,
θ
t
)
=
∑
i
=
1
N
∑
z
i
=
C
1
,
.
.
.
,
C
K
[
l
o
g
(
p
(
x
i
,
z
i
∣
θ
)
)
p
(
z
i
∣
x
i
,
θ
t
)
]
Q(\theta,\theta^{t})=\sum_{i=1}^{N}\sum_{z_{i}=C_{1},...,C_{K}}[log(p(x_{i},z_{i}|\theta))p(z_{i}|x_{i},\theta^{t})]
Q(θ,θt)=i=1∑Nzi=C1,...,CK∑[log(p(xi,zi∣θ))p(zi∣xi,θt)]
=
∑
i
=
1
N
∑
z
i
=
C
1
,
.
.
.
,
C
K
[
l
o
g
(
p
z
i
N
(
x
i
∣
μ
z
i
,
Σ
z
i
)
)
]
p
z
i
t
N
(
x
i
∣
μ
z
i
t
,
Σ
z
i
t
)
∑
k
=
1
K
p
k
t
N
(
x
i
∣
μ
k
t
,
Σ
k
t
)
=\sum_{i=1}^{N}\sum_{z_{i}=C_{1},...,C_{K}}[log(p_{z_{i}}N(x_{i}|\mu_{z_{i}},\Sigma_{z_{i}}))]\frac{p_{z_{i}}^{t}N(x_{i}|\mu_{z_{i}}^{t},\Sigma_{z_{i}}^{t})}{\sum_{k=1}^{K}p_{k}^{t}N(x_{i}|\mu_{k}^{t},\Sigma_{k}^{t})}
=i=1∑Nzi=C1,...,CK∑[log(pziN(xi∣μzi,Σzi))]∑k=1KpktN(xi∣μkt,Σkt)pzitN(xi∣μzit,Σzit)
这就是在第
t
t
t轮向第
t
+
1
t+1
t+1轮迭代的E步，所有带上标
t
t
t的参数都是已知的，都是上一步的迭代结果。下一步需要执行M步：
θ
t
+
1
=
a
r
g
m
a
x
θ
Q
(
θ
,
θ
t
)
\theta^{t+1}=argmax_{\theta}Q(\theta,\theta^{t})
θt+1=argmaxθQ(θ,θt)
参数的估计结果如下：
p
k
t
+
1
=
1
N
∑
i
=
1
N
p
(
z
i
=
C
k
∣
x
i
,
θ
t
)
p_{k}^{t+1}=\frac{1}{N}\sum_{i=1}^{N}p(z_{i}=C_{k}|x_{i},\theta^{t})
pkt+1=N1i=1∑Np(zi=Ck∣xi,θt)
μ
k
t
+
1
=
∑
i
=
1
N
p
(
z
i
=
C
k
∣
x
i
,
θ
t
)
x
i
∑
i
=
1
N
p
(
z
i
=
C
k
∣
x
i
,
θ
t
)
\mu_{k}^{t+1}=\frac{\sum_{i=1}^{N}p(z_{i}=C_{k}|x_{i},\theta^{t})x_{i}}{\sum_{i=1}^{N}p(z_{i}=C_{k}|x_{i},\theta^{t})}
μkt+1=∑i=1Np(zi=Ck∣xi,θt)∑i=1Np(zi=Ck∣xi,θt)xi
Σ
k
t
+
1
=
∑
i
=
1
N
p
(
z
i
=
C
k
∣
x
i
,
θ
t
)
(
x
i
−
μ
k
t
+
1
)
(
x
i
−
μ
k
t
+
1
)
T
∑
i
=
1
N
p
(
z
i
=
C
k
∣
x
i
,
θ
t
)
\Sigma_{k}^{t+1}=\frac{\sum_{i=1}^{N}p(z_{i}=C_{k}|x_{i},\theta^{t})(x_{i}-\mu_{k}^{t+1})(x_{i}-\mu_{k}^{t+1})^{T}}{\sum_{i=1}^{N}p(z_{i}=C_{k}|x_{i},\theta^{t})}
Σkt+1=∑i=1Np(zi=Ck∣xi,θt)∑i=1Np(zi=Ck∣xi,θt)(xi−μkt+1)(xi−μkt+1)T
高斯混合模型的应用场景
当得到高斯混合模型的参数后，将可以判断样本属于哪一类（哪一个高斯分布），分别计算样本属于 K K K个高斯分布的概率，选择概率值最高的高斯分布作为其类别。
样本
x
i
x_{i}
xi属于第
k
k
k个高斯分布的概率：
p
(
z
i
=
C
k
∣
x
i
)
=
p
(
z
i
=
C
k
)
p
(
x
i
∣
z
i
=
C
k
)
p
(
x
i
)
∝
p
(
z
i
=
C
k
)
p
(
x
i
∣
z
i
=
C
k
)
=
p
k
N
(
x
i
∣
μ
k
,
Σ
k
)
p(z_{i}=C_{k}|x_{i})=\frac{p(z_{i}=C_{k})p(x_{i}|z_{i}=C_{k})}{p(x_{i})}\propto p(z_{i}=C_{k})p(x_{i}|z_{i}=C_{k})=p_{k}N(x_{i}|\mu_{k},\Sigma_{k})
p(zi=Ck∣xi)=p(xi)p(zi=Ck)p(xi∣zi=Ck)∝p(zi=Ck)p(xi∣zi=Ck)=pkN(xi∣μk,Σk)
样本
x
i
x_{i}
xi对于每个高斯分布都能计算出一个上面形式的概率值，取最大的那一个作为样本属于的高斯分布，即得到样本的类别。
以下是实验演示：
%matplotlib inline
import matplotlib.pyplot as plt
from sklearn.mixture import GaussianMixture
"""
生成多个各向同性的高斯分布：
make_blobs(
n_samples=100, # 所有簇的样本数
n_features=2, # 每个样本的特征数
centers=None, # 簇的数量
cluster_std=1.0, # 每个簇的协方差矩阵元素值
center_box=(-10.0, 10.0), # 随机生成簇时每个簇的边界框。
shuffle=True,
random_state=None)
各向同性高斯分布（球形高斯分布）：各个方向方差都一样的多维高斯分布
协方差矩阵Sigma=数值sigma*单位阵I
"""
from sklearn.datasets.samples_generator import make_blobs
# 产生实验数据：
# X：(1000,2)
# y_true：(1000,)取值范围0,1,2,3
X, y_true = make_blobs(n_samples=1000, centers=4)
# sharey='row'每个子图的行将共享同一个y坐标
# ax 为坐标系
fig, ax = plt.subplots(1,2,sharey='row')
# s：刻度单位；alpha：透明度
ax[0].scatter(X[:, 0], X[:, 1], s=5, alpha=0.5)
#高斯混合模型拟合样本
gmm = GaussianMixture(n_components=4)
gmm.fit(X)
#打印GMM三组参数：权重、均值、协方差矩阵
print(gmm.weights_)
print(gmm.means_)
print(gmm.covariances_)
#打印前10个样本属于每个高斯分布的概率（保留5位小数）
print(gmm.predict_proba(X)[:10].round(5))
#通过GMM模型推测每个样本所属的类别
labels = gmm.predict(X)
print(labels)
#不同的类别标记为不同的颜色
#cmap='viridis':黄到蓝
ax[1].scatter(X[:, 0], X[:, 1], s=5, alpha=0.5, c=labels, cmap='viridis')
plt.show()
结果为：
[0.25156792 0.25020204 0.25072685 0.24750319]
[[-6.2715644 -3.75712384]
[ 0.19631038 2.40819993]
[ 2.10422918 -7.32369779]
[-2.06228862 -3.91610721]]
[[[ 1.02207524 0.01495611]
[ 0.01495611 1.0354908 ]]
[[ 0.98060509 -0.01132009]
[-0.01132009 1.00077596]]
[[ 0.90053266 0.11036538]
[ 0.11036538 1.0409055 ]]
[[ 0.90831192 -0.01799792]
[-0.01799792 0.85407048]]]
[[0.0000e+00 0.0000e+00 1.0000e+00 0.0000e+00]
[1.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00]
[0.0000e+00 1.0000e+00 0.0000e+00 0.0000e+00]
[0.0000e+00 0.0000e+00 1.0000e+00 0.0000e+00]
[6.1800e-03 0.0000e+00 0.0000e+00 9.9382e-01]
[4.2300e-03 0.0000e+00 0.0000e+00 9.9577e-01]
[0.0000e+00 0.0000e+00 1.0000e+00 0.0000e+00]
[0.0000e+00 0.0000e+00 1.0000e+00 0.0000e+00]
[1.0000e-05 0.0000e+00 2.0000e-05 9.9997e-01]
[2.0000e-05 0.0000e+00 0.0000e+00 9.9998e-01]]
[2 0 1 2 3 3 2 2 3 3 1 2 1 0 2 3 0 0 2 2 2 3 0 1 0 2 1 0 3 1 2 0 2 0 1 0 1
1 0 0 2 1 1 3 2 3 0 2 0 0 3 3 3 3 3 1 2 1 1 1 1 0 1 1 2 0 2 2 3 1 3 3 3 0
2 3 3 1 1 1 1 2 0 1 0 1 3 2 1 0 0 3 2 0 3 0 1 0 2 0 1 1 1 0 3 0 1 0 2 2 3
3 3 2 2 3 0 3 1 3 0 0 3 2 0 3 2 3 3 3 2 0 1 1 2 3 3 0 0 0 2 1 2 1 2 0 0 0
2 3 3 2 1 0 1 0 1 2 1 1 2 2 0 2 3 1 1 1 2 2 3 2 3 2 2 1 2 3 1 0 0 3 3 2 0
0 1 1 0 2 3 3 2 2 3 3 2 2 3 1 0 0 1 3 2 0 1 2 1 1 2 1 3 3 2 1 3 2 2 1 3 1
0 3 1 3 0 0 2 0 1 2 0 3 3 3 1 3 3 0 0 1 0 0 3 2 2 2 2 3 1 3 1 1 0 2 0 0 3
0 2 1 2 1 1 1 0 2 1 1 1 2 0 2 0 0 3 0 2 1 2 2 0 2 2 3 1 1 0 3 2 2 0 3 1 3
0 3 0 0 3 1 3 3 0 2 1 0 3 3 1 2 2 1 1 3 1 2 2 0 3 2 0 3 0 0 0 1 2 3 1 3 2
3 3 3 3 0 2 0 1 3 3 2 2 1 1 3 3 1 0 2 2 3 3 3 3 0 2 0 1 1 2 3 2 0 0 1 1 2
1 3 1 1 3 3 2 3 3 1 1 1 1 0 0 1 2 3 2 3 3 0 0 3 1 0 0 3 0 1 1 3 2 0 2 2 0
3 3 3 3 3 2 3 2 3 3 2 1 0 3 0 1 0 1 1 3 0 1 2 0 3 0 0 2 1 1 0 1 0 2 3 3 2
3 3 3 0 3 1 0 0 1 3 2 3 1 1 2 0 0 1 1 1 1 0 2 2 1 3 2 2 2 1 1 3 1 3 2 1 0
1 0 1 0 3 1 2 2 0 2 3 1 2 3 1 3 1 1 2 0 2 1 3 0 2 0 2 2 0 3 2 1 1 0 1 1 1
0 0 2 3 3 3 2 0 0 0 3 3 0 0 2 3 1 0 1 3 0 2 1 1 2 3 1 2 2 0 0 0 2 0 3 0 0
2 2 2 1 3 2 0 0 1 3 3 1 1 1 2 2 0 0 2 1 3 1 1 2 1 0 2 1 1 2 3 1 0 1 1 0 1
0 1 1 3 2 1 1 0 0 3 0 2 3 1 2 2 2 0 2 0 3 1 0 3 0 1 0 2 3 3 0 2 3 3 0 2 0
3 3 0 1 2 2 2 1 2 0 2 1 0 3 3 2 1 1 2 1 2 3 0 3 1 3 3 3 0 1 0 0 0 0 2 0 3
0 1 2 1 1 1 2 0 3 0 2 3 1 3 3 1 2 2 1 0 2 2 3 3 3 1 0 2 2 1 0 1 1 2 0 0 0
2 2 2 2 1 0 1 2 3 2 1 1 1 0 3 1 2 3 1 1 3 2 2 2 3 0 3 1 2 2 2 0 3 2 1 3 3
2 0 3 2 0 1 0 0 3 0 2 0 0 2 2 1 0 1 2 2 1 3 0 0 3 3 0 0 3 2 2 0 3 3 1 0 0
1 0 0 3 3 3 2 3 3 1 0 3 2 3 1 0 3 3 1 0 3 3 2 0 1 1 2 0 1 3 1 2 3 1 1 2 1
0 0 2 3 1 0 1 2 3 3 1 0 0 2 0 0 3 3 0 1 3 0 2 1 0 2 1 0 2 2 1 0 2 0 2 0 2
2 0 0 2 2 3 3 0 2 1 2 1 3 0 1 2 0 2 2 3 3 3 1 3 2 3 3 2 0 3 1 2 2 0 2 2 3
3 2 0 0 2 2 2 1 0 2 0 2 2 1 0 0 3 0 0 0 3 2 3 2 1 2 2 2 1 1 0 0 0 1 2 2 3
3 0 0 1 0 1 3 2 1 2 3 2 3 0 3 1 1 3 3 0 3 1 3 2 3 0 3 1 2 0 1 1 3 1 1 2 1
1 3 0 1 1 1 3 1 2 0 1 1 0 1 1 0 2 3 1 0 3 3 3 1 1 3 0 3 2 2 2 3 0 0 0 2 0
3]
扫一扫
1540
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
