李一鸣
1160300625
2018 年 9 月 28 日
生成随机序列 $\mathrm{\mathbf{X}} = (X_1, X_2, ..., X_n)$，其中每个 $X_i$ 服从 $[-\frac{a}{2}, \frac{a}{2}]$ 的均匀分布。
生成随机序列 $\mathrm{\mathbf{Y}} = (Y_1, Y_2, ..., Y_n)$，其中每个 $Y_i$ 服从 $[-\frac{a}{2}, \frac{a}{2}]$ 的均匀分布。
蒙特卡罗投点法：
在边长为 $a$ 的正方形内随机投点，设该点落入此正方形的内切圆中的概率为 $P_{circle}$，则：
$$
\begin{aligned}
P_{circle} &= \frac{S_{circle}}{S_{square}} \\
&= \frac{\pi (\frac{a}{2})^2}{a^2} \\
&= \frac{\pi}{4}
\end{aligned}
\tag{1}
$$
假定总共生成了 $n$ 个数据，其中有 $m$ 个在圆内，则：
$$
f_{circle} = \frac{m}{n}
\tag{2}
$$
对于任一点 $(X_i, Y_i)$，如果满足：
$$
X_i^2 + Y_i^2 \le (\frac{a}{2})^2 = \frac{a^2}{4}
\tag{3}
$$
则其在圆内，计入 $m$ 中。
以频率估计概率，我们有：
$$
\pi = \frac{4m}{n}
\tag{4}
$$
实验结果：
在实验中取 $a = 1$（其实 $a$ 取多少都没有关系，精确度只与样本数相关），取 $n = 1\ 000, 10\ 000, 100\ 000, 1\ 000\ 000$ 分别进行实验。
生成均值为 $\mu = 10$、方差为 $\sigma = 5$ 的正态分布，并画出均值和方差随样本数增加而变化的图。
设样本为 $\mathrm{\mathbf{X}}$，总样本数为 $N$，记前 $n$ 个样本数据的均值、方差分别为 $E_n$ 和 $D_n$，则得到均值、方差矩阵：
$$
\mathrm{\mathbf{E}} = (E_1, E_2, ..., E_N) \\
\mathrm{\mathbf{D}} = (D_1, D_2, ..., D_N)
\tag{5}
$$
样本个数矩阵：
$$
\mathrm{\mathbf{N}} = (1, 2, ..., N)
\tag{6}
$$
我们只需要作出 $(\mathrm{\mathbf{N}}, \mathrm{\mathbf{E}})$ 和 $(\mathrm{\mathbf{N}}, \mathrm{\mathbf{D}})$ 的图像即可。
注意事项：
本来我们可以直接根据 $\mathrm{\mathbf{X}}$ 中的前 $n$ 项直接计算 $E_n$ 和 $D_n$，但在实际问题中数据量往往过大，我们一般不会保存之前的数据，因此我们通常采用递推的计算方式。
-
均值递推公式
$$
\begin{cases}
E_{n - 1} = \frac{1}{n - 1}(X_1 + X_2 + ... + X_{n-1}) \ \
E_{n} = \frac{1}{n}(X_1 + X_2 + ... + X_{n - 1} + X_{n}) \
\end{cases}
$$
解得：
$$
\begin{aligned}
E_n &= \frac{(n - 1)E_{n - 1} + X_n}{n} \
&= \frac{nE_{n - 1} - E_{n - 1} + X_n}{n} \
&= E_{n - 1} + \frac{X_n - E_{n - 1}}{n}
\tag{ps.1}
\end{aligned}
$$
-
方差递推公式
$$
\begin{cases}
D_{n - 1} = \frac{1}{n - 1}\sum_{i = 1}^{n - 1}(X_i - E_{n - 1})^2 \ \
D_{n} = \frac{1}{n}\sum_{i = 1}^{n}(X_i - E_{n})^2 \
\end{cases}
$$
联立式 (ps.1)，得：
$$
\begin{aligned}
D_{n} &= \frac{1}{n}\sum_{i = 1}^n(X_i - E_{n - 1} - \frac{X_n - E_{n - 1}}{n})^2 \
&= \frac{1}{n}\sum_{i = 1}^n[(X_i - E_{n - 1})^2 + (\frac{X_n - E_{n - 1}}{n})^2 - 2(X_i - E_{n - 1})(\frac{X_n - E_{n - 1}}{n})] \
&= (\frac{X_n - E_{n - 1}}{n})^2 + \frac{1}{n}\sum_{i = 1}^n[(X_i - E_{n - 1})^2 - 2(X_i - E_{n - 1})(\frac{X_n - E_{n - 1}}{n})] \
&= (\frac{X_n - E_{n - 1}}{n})^2 + \frac{1}{n}[(X_n - E_{n - 1})^2 - 2(X_n - E_{n - 1})(\frac{X_n - E_{n - 1}}{n})] \
& \quad \quad \quad \quad \quad \quad \quad \quad + \frac{1}{n}\sum_{i = 1}^{n - 1}[(X_i - E_{n - 1})^2 - 2(X_i - E_{n - 1})(\frac{X_n - E_{n - 1}}{n})] \
&= \frac{n - 1}{n^2}(X_n - E_{n - 1})^2 + \frac{1}{n}\sum_{i = 1}^{n - 1}(X_i - E_{n - 1})^2 - 2(\frac{X_n - E_{n - 1}}{n})\sum_{i = 1}^{n - 1}(X_i - E_{n - 1}) \
&= \frac{n - 1}{n^2}(X_n - E_{n - 1})^2 + \frac{n - 1}{n}\frac{1}{n - 1}\sum_{i = 1}^{n - 1}(X_i - E_{n - 1})^2) \
&= \frac{n - 1}{n^2}(X_n - D_{n - 1})^2 + \frac{n - 1}{n}D_{n - 1}
\tag{ps.2}
\end{aligned}
$$
实验结果：
在实验中取 $a = 1$，取 $N = 1\ 000, 10\ 000$ 分别进行实验。
敌军坦克分队到达我方阵地规律服从泊松分布，平均每分钟到达 $\lambda$ 辆。
泊松分布的期望值是 $\lambda$，也就是说在一分钟之内，到达的坦克数量 $T$ 的分布列为：
$$
P(T = k) = \frac{\lambda ^ k e^{-\lambda}}{k!}
\tag{7}
$$
我们可以生成 $N$ 组数据 $\mathrm{\mathbf{T_1}}, \mathrm{\mathbf{T_2}}, ..., \mathrm{\mathbf{T_N}}$，分别用它们的均值 $E(\mathrm{\mathbf{T_1}}), E(\mathrm{\mathbf{T_2}}),..., E(\mathrm{\mathbf{T_N}})$ 表示第 $1, 2, ..., N$ 分钟内到达的坦克数量存入 $\mathrm{\mathbf{A}} = (A_1, A_2, ..., A_N) 中$，则在 $N$ 分钟内坦克到达总数量 $A_{total}$ 满足：
$$
A_{total} = \sum_{n = 1}^N A_n
\tag{8}
$$
实验结果：
取 $\lambda = 4, N = 3$，对样本数进行改变，得到如下实验结果：
每辆敌军坦克到达的时刻服从期望为 $\frac{1}{\lambda}$ 的指数分布，也就是说坦克到达的时间 $S$ 的分布函数为：
$$
F_S(x) = e^{-\lambda x}
\tag{9}
$$
我们可以生成 $M$ 组数据 $\mathrm{\mathbf{S_1}}, \mathrm{\mathbf{S_2}}, ..., \mathrm{\mathbf{S_M}}$，分别用它们的均值 $E(\mathrm{\mathbf{S_1}}), E(\mathrm{\mathbf{S_2}}),..., E(\mathrm{\mathbf{S_M}})$ 表示第 $1, 2, ..., M$ 辆坦克到达所需的时间，将其存入 $\mathrm{\mathbf{B}} = (B_1, B_2, ..., B_M) 中$，则在 $N$ 分钟内每辆敌军坦克到达时间为：
$$
\mathrm{\mathbf{B'}} = (B'_i = \sum_1^j B_j | \quad j \in [1, M] \quad where \quad B'_i < N)
\tag{10}
$$
实验结果：
取 $\lambda = 4, N = 3$，对样本组数 $M$ 进行改变，得到如下实验结果：
考虑水果聚类问题，水果的属性 $\mathrm{\mathbf{X}}$ 满足高斯分布，其均值向量、协方差矩阵分别为 $\mathrm{\mathbf{\mu}}, \mathrm{\mathbf{\Sigma}}$，将其概率密度记为 $p(\mathrm{\mathbf{X}}|\mathrm{\mathbf{\mu}}, \mathrm{\mathbf{\Sigma}})$。
定义高斯混合分布：
$$
\begin{aligned}
p_M(\mathrm{\mathbf{X}}) &= \sum_{i = 1}^k \alpha_i p(\mathrm{\mathbf{X}}|\mathrm{\mathbf{\mu _i}}, \mathrm{\mathbf{\Sigma _i}}) \\
p(\mathrm{\mathbf{X}}|\mathrm{\mathbf{\mu _i}}, \mathrm{\mathbf{\Sigma _i}}) &= \frac{exp{-\frac{1}{2}(\mathrm{\mathbf{X}} - \mathrm{\mathbf{\mu_i}})^T\mathrm{\mathbf{\Sigma_i^{-1}}}(\mathrm{\mathbf{X}} - \mathrm{\mathbf{\mu_i}})}}{(2\pi) ^ {D/2}|\Sigma_i|^{\frac{1}{2}}}
\end{aligned}
\tag{11}
$$
该分布由 $k$ 个混合分布组成，每个混合成分对应一个高斯分布，其中 $\mu _i, \Sigma _i$ 是第 $i$ 个高斯混合成分的参数，$\mathrm{\mathbf{\alpha _i}}$ 为选择第 $i$ 个混合成分的概率，满足：
$$
\alpha i > 0, \sum{i = 1}^k \alpha_i = 1
\tag{12}
$$
记样本 $X_j$ 属于第 $i$ 个高斯成分的后验概率为 $y_{ji}$，有：
$$
\begin{aligned}
y_{ji} &= \frac{\alpha _i p(x_j|\mu _i, \Sigma _i)}{p_M(x_j)} \
&= \frac{\alpha _i p(x_j|\mu _i, \Sigma i)}{\sum{l = 1}^k \alpha _l p(x_j|\mu _l, \Sigma_l)}
\tag{13}
\end{aligned}
$$
为了得到混合分布的各个组成部分的分布参数，我们需要利用 EM 算法 (Expectation–maximization algorithm) 不断迭代来获取 $k$ 个类的均值和方差参数。
E 步：
根据当前参数计算样本后验概率 $\mathrm{\mathbf{Y}} = (y_{ji})_{ji}$
M 步：
根据后验概率更新模型参数 ${\alpha _i, \mu _i, \Sigma _i | 1 \le i \le k}$，新的参数与后验概率应该满足下面的关系：
$$
\begin{aligned}
\alpha_i' &= \frac{\sum {j = 1}^N y{ji}}{N} \ \
\mathrm{\mathbf{\mu_i'}} &= \frac{\sum {j = 1}^N y{ji}x_j}{\sum {j = 1}^N y{ji}} \ \
\mathrm{\mathbf{\Sigma_i'}} &= \frac{\Sigma_{j = 1}^{N}y_{ji}(x_j - \mu_i')(x_j - \mu_i')^T}{\sum {j = 1}^{N}y{ji}}
\tag{14}
\end{aligned}
$$
不断重复 E、M 两步直到收敛。
实验结果：
现有水果数据 $\mathrm{\mathbf{S}}$：
$$
\mathrm{\mathbf{S}} = (\mathrm{\mathbf{S_1}}, \mathrm{\mathbf{S_2}}, ..., \mathrm{\mathbf{S_{N}}})
\tag{15}
$$
其中 N = 30，$S_i$ 为二维列向量，包含密度、含糖率两个属性，我们随机初始化一组参数：
$$
\begin{aligned}
\alpha_1 &= \alpha_2 = \alpha_3 = \frac{1}{3} \\
\mathrm{\mathbf{\mu_1}} &= \mathrm{\mathbf{S_6}}, \mathrm{\mathbf{\mu_2}} = \mathrm{\mathbf{S_{22}}}, \mathrm{\mathbf{\mu_3}} = \mathrm{\mathbf{S_{27}}} \\
\mathrm{\mathbf{\Sigma_1}} &= \mathrm{\mathbf{\Sigma_2}} = \mathrm{\mathbf{\Sigma_3}} = \begin{pmatrix}
0.1 & 0.0 \\
0.0 & 0.1 \\
\end{pmatrix}
\tag{16}
\end{aligned}
$$
令迭代次数 $I = 50$，将得到的结果以散点图绘制如下：
详细计算过程参见 em-50.txt。
- [Monte Carlo method | Wikipedia]
- [Expectation–maximization algorithm | Wikipedia]
