现代机器学习（特别是大规模数据场景）中，往往用随机模拟和统计推断配合进行处理：
用统计推断建立模型 ：比如用EM算法算出一个高斯混合模型（GMM），找到了数据的分布规律（参数）。用随机模拟验证或使用模型 ：得出规律后，我们需要生成几百万个模拟数据来测试这个模型好不好用；或者在实际业务中，利用这个算出来的分布规律，随机模拟用户未来的点击行为。分布式加速 ：无论是模拟生成上亿个随机数，还是对海量数据进行EM算法迭代，单机都算不动，必须依靠 Spark 把数据和计算拆分到多台机器上并行处理。随机模拟是已知公式，让Spark集群算出上亿级别的模拟样本；
统计推断是已知这上亿个样本，让Spark集群算出背后的公式参数。
相关的方法主要如下：
三种随机数生成法 :
分布函数的逆写得出来 → 逆分布法;
写不出来但能找到试投密度 g g g c c c 逆分布法公式 :
X = F − 1 ( U ) X = F^{-1}(U) X = F − 1 ( U ) U ∼ U ( 0 , 1 ) U\sim U(0,1) U ∼ U ( 0 , 1 ) X = − 1 λ ln U X=-\tfrac{1}{\lambda}\ln U X = − λ 1 ln U 拒绝接受法 :
接受条件 U ≤ f ( Y ) c g ( Y ) U \le \dfrac{f(Y)}{c\,g(Y)} U ≤ c g ( Y ) f ( Y ) = 1 c =\dfrac{1}{c} = c 1 c c c EM 四步 :
选初值 → E 步算 Q 函数(隐变量的条件期望)→ M 步极大化 Q 更新参数 → 迭代到收敛。遗传模型迭代式 :
z ( i ) = 125 ⋅ θ ( i ) 2 + θ ( i ) z^{(i)} = 125\cdot\dfrac{\theta^{(i)}}{2+\theta^{(i)}} z ( i ) = 125 ⋅ 2 + θ ( i ) θ ( i ) θ ( i + 1 ) = z ( i ) + 34 z ( i ) + 72 \theta^{(i+1)} = \dfrac{z^{(i)}+34}{z^{(i)}+72} θ ( i + 1 ) = z ( i ) + 72 z ( i ) + 34 GMM 三个更新式 :
μ k \mu_k μ k σ k 2 \sigma_k^2 σ k 2 α k \alpha_k α k γ ^ j k \hat\gamma_{jk} γ ^ j k 分布式 EM :
充分统计量可加,各节点局部求和、全局汇总更新参数计算机只会生成均匀分布的随机数，而我们生成符合对应分布的随机数需要进行一个映射关系的处理，逆分布法是最直接的一种随机数生成方法。
当随机变量所服从的分布函数已知,并且这个分布函数的逆函数有显式解的时候,就可以直接用它来产生服从该分布的随机数。一旦分布函数的逆写不出来,这个方法就用不了,得改用后面的拒绝接受法。
原理 :
设随机变量 X X X F ( x ) F(x) F ( x ) U U U [ 0 , 1 ] [0,1] [ 0 , 1 ] X = F − 1 ( U ) X = F^{-1}(U) X = F − 1 ( U ) F F F
算法步骤 :
产生一个服从均匀分布 U ∼ U ( 0 , 1 ) U \sim U(0,1) U ∼ U ( 0 , 1 ) 令 X = F − 1 ( U ) X = F^{-1}(U) X = F − 1 ( U ) X X X F F F 要产生 n n n Exp ( λ ) \text{Exp}(\lambda) Exp ( λ )
f ( x ) = λ e − λ x , F ( x ) = 1 − e − λ x ( x ≥ 0 ) f(x) = \lambda e^{-\lambda x}, \quad F(x) = 1 - e^{-\lambda x} \quad (x \ge 0) f ( x ) = λ e − λ x , F ( x ) = 1 − e − λ x ( x ≥ 0 )
令 U = F ( X ) = 1 − e − λ X U = F(X) = 1 - e^{-\lambda X} U = F ( X ) = 1 − e − λ X
X = − 1 λ ln ( 1 − U ) X = -\frac{1}{\lambda}\ln(1 - U) X = − λ 1 ln ( 1 − U )
1 − U 1-U 1 − U U U U U ( 0 , 1 ) U(0,1) U ( 0 , 1 ) X = − 1 λ ln U X = -\dfrac{1}{\lambda}\ln U X = − λ 1 ln U
要产生 n n n
f ( x ) = 1 2 e − ∣ x ∣ f(x) = \frac{1}{2} e^{-|x|} f ( x ) = 2 1 e − ∣ x ∣
对它积分,可得累积分布函数为分段形式:
F ( x ) = { 1 2 e x , x < 0 1 − 1 2 e − x , x ≥ 0 F(x) = \begin{cases} \dfrac{1}{2}e^{x}, & x < 0 \\[2mm] 1 - \dfrac{1}{2}e^{-x}, & x \ge 0 \end{cases} F ( x ) = ⎩ ⎨ ⎧ 2 1 e x , 1 − 2 1 e − x , x < 0 x ≥ 0
令 U = F ( X ) U = F(X) U = F ( X )
X = { ln ( 2 U ) , 0 < U < 1 2 − ln ( 2 ( 1 − U ) ) , 1 2 ≤ U < 1 X = \begin{cases} \ln(2U), & 0 < U < \tfrac{1}{2} \\[2mm] -\ln\big(2(1-U)\big), & \tfrac{1}{2} \le U < 1 \end{cases} X = ⎩ ⎨ ⎧ ln ( 2 U ) , − ln ( 2 ( 1 − U ) ) , 0 < U < 2 1 2 1 ≤ U < 1
按这个式子把均匀随机数代进去,就得到 Laplace 分布的随机数。
当要生成的随机数的分布没有分布函数,或者分布函数的逆没有显式解 时,逆分布法就失效了,这时候用拒绝接受法。
适用条件 :能够找到另一个概率密度函数 g ( x ) g(x) g ( x ) f ( x ) f(x) f ( x ) c c c
f ( x ) ≤ c g ( x ) f(x) \le c\,g(x) f ( x ) ≤ c g ( x )
这里 g g g 试投密度(建议分布) ,属于容易直接抽样的分布(比如均匀分布、正态分布)。
算法步骤 :
产生一个密度函数为 g g g Y Y Y U ∼ U ( 0 , 1 ) U \sim U(0,1) U ∼ U ( 0 , 1 ) 如果满足 U ≤ f ( Y ) c g ( Y ) U \le \frac{f(Y)}{c\,g(Y)} U ≤ c g ( Y ) f ( Y ) X = Y X = Y X = Y 推导思路 :令接受条件成立的那个事件,可以证明"被接受的 Y Y Y f f f Y Y Y ∫ f \int f ∫ f f f f
P ( 接受 ∣ Y = y ) = f ( y ) c g ( y ) , P ( 接受 ) = 1 c P(\text{接受} \mid Y=y) = \frac{f(y)}{c\,g(y)}, \qquad P(\text{接受}) = \frac{1}{c} P ( 接受 ∣ Y = y ) = c g ( y ) f ( y ) , P ( 接受 ) = c 1
拒绝接受算法的效率由接受率 决定,而接受率正好是 1 c \dfrac{1}{c} c 1 c c c c c c
舍选法Ⅰ是拒绝接受法在一种简单情形下的特例
适用条件 :要生成的随机数的分布取值于一个有限区间 [ a , b ] [a,b] [ a , b ] f ( x ) f(x) f ( x ) M M M x x x f ( x ) ≤ M f(x) \le M f ( x ) ≤ M
算法步骤 :
在区间上生成 X 1 ∼ U ( a , b ) X_1 \sim U(a,b) X 1 ∼ U ( a , b ) X 2 ∼ U ( 0 , M ) X_2 \sim U(0, M) X 2 ∼ U ( 0 , M ) 如果 X 2 ≤ f ( X 1 ) X_2 \le f(X_1) X 2 ≤ f ( X 1 ) X = X 1 X = X_1 X = X 1 这实际上就是在矩形框 [ a , b ] × [ 0 , M ] [a,b] \times [0,M] [ a , b ] × [ 0 , M ] f f f
产生随机数 X 1 X_1 X 1 U ∼ U ( 0 , 1 ) U \sim U(0,1) U ∼ U ( 0 , 1 ) 如果 U ≤ f ( X 1 ) M U \le \dfrac{f(X_1)}{M} U ≤ M f ( X 1 ) X = X 1 X = X_1 X = X 1 两种写法本质是同一件事,只是把 X 2 = M ⋅ U X_2 = M\cdot U X 2 = M ⋅ U
EM(Expectation-Maximization) 算法用来求含 隐变量 的概率模型的参数极大似然估计。当似然函数因为存在隐变量而难以直接极大化时,EM 通过"E 步补隐变量的期望、M 步再极大化"两步交替迭代,逐步逼近最优解。
算法步骤 :
选初值 :选择参数的初值 θ ( 0 ) \theta^{(0)} θ ( 0 ) E 步 :记 θ ( i ) \theta^{(i)} θ ( i ) i i i i + 1 i+1 i + 1 Q ( θ , θ ( i ) ) = E Z ∣ Y , θ ( i ) [ log P ( Y , Z ∣ θ ) ] = ∑ Z log P ( Y , Z ∣ θ ) P ( Z ∣ Y , θ ( i ) ) Q(\theta, \theta^{(i)}) = E_{Z\mid Y,\theta^{(i)}}\big[\log P(Y,Z \mid \theta)\big] = \sum_{Z} \log P(Y,Z\mid\theta)\, P(Z\mid Y,\theta^{(i)}) Q ( θ , θ ( i ) ) = E Z ∣ Y , θ ( i ) [ log P ( Y , Z ∣ θ ) ] = ∑ Z log P ( Y , Z ∣ θ ) P ( Z ∣ Y , θ ( i ) ) P ( Z ∣ Y , θ ( i ) ) P(Z\mid Y,\theta^{(i)}) P ( Z ∣ Y , θ ( i ) ) Y Y Y θ ( i ) \theta^{(i)} θ ( i ) Z Z Z M 步 :求使 Q ( θ , θ ( i ) ) Q(\theta,\theta^{(i)}) Q ( θ , θ ( i ) ) θ \theta θ i + 1 i+1 i + 1 θ ( i + 1 ) = arg max θ Q ( θ , θ ( i ) ) \theta^{(i+1)} = \arg\max_{\theta} Q(\theta, \theta^{(i)}) θ ( i + 1 ) = arg max θ Q ( θ , θ ( i ) ) 迭代 :重复第 2、3 步,直至收敛(参数或对数似然的变化足够小)。设一次实验可能出现四个结果,发生概率分别为
1 2 + θ 4 , 1 − θ 4 , 1 − θ 4 , θ 4 , 0 ≤ θ ≤ 1 \frac{1}{2} + \frac{\theta}{4},\quad \frac{1-\theta}{4},\quad \frac{1-\theta}{4},\quad \frac{\theta}{4}, \qquad 0 \le \theta \le 1 2 1 + 4 θ , 4 1 − θ , 4 1 − θ , 4 θ , 0 ≤ θ ≤ 1
现在进行了 197 次试验,四种结果的发生次数分别为 125、18、20、34 。要求出分布中参数 θ \theta θ
引入隐变量 :第一类结果的概率 1 2 + θ 4 \tfrac{1}{2} + \tfrac{\theta}{4} 2 1 + 4 θ 1 2 \tfrac{1}{2} 2 1 θ 4 \tfrac{\theta}{4} 4 θ y 1 = 125 − z y_1 = 125 - z y 1 = 125 − z z z z 1 2 \tfrac{1}{2} 2 1 θ 4 \tfrac{\theta}{4} 4 θ z z z
完全似然函数 可以写成(带隐变量 z z z
L ( θ ) ∝ ( 1 2 ) 125 − z ( θ 4 ) z ( 1 − θ 4 ) 18 ( 1 − θ 4 ) 20 ( θ 4 ) 34 L(\theta) \propto \left(\frac{1}{2}\right)^{125-z}\left(\frac{\theta}{4}\right)^{z}\left(\frac{1-\theta}{4}\right)^{18}\left(\frac{1-\theta}{4}\right)^{20}\left(\frac{\theta}{4}\right)^{34} L ( θ ) ∝ ( 2 1 ) 125 − z ( 4 θ ) z ( 4 1 − θ ) 18 ( 4 1 − θ ) 20 ( 4 θ ) 34
对数似然 (去掉与 θ \theta θ
log L ( θ ) ∝ ( z + 34 ) log θ + ( 18 + 20 ) log ( 1 − θ ) \log L(\theta) \propto (z + 34)\log\theta + (18+20)\log(1-\theta) log L ( θ ) ∝ ( z + 34 ) log θ + ( 18 + 20 ) log ( 1 − θ )
E 步 :在当前估计 θ ( i ) \theta^{(i)} θ ( i ) z z z
z ( i ) = E [ z ∣ θ ( i ) ] = 125 ⋅ θ ( i ) 4 1 2 + θ ( i ) 4 = 125 ⋅ θ ( i ) 2 + θ ( i ) z^{(i)} = E\big[z \mid \theta^{(i)}\big] = 125 \cdot \frac{\tfrac{\theta^{(i)}}{4}}{\tfrac{1}{2} + \tfrac{\theta^{(i)}}{4}} = 125 \cdot \frac{\theta^{(i)}}{2 + \theta^{(i)}} z ( i ) = E [ z ∣ θ ( i ) ] = 125 ⋅ 2 1 + 4 θ ( i ) 4 θ ( i ) = 125 ⋅ 2 + θ ( i ) θ ( i )
M 步 :对 Q Q Q z z z z ( i ) z^{(i)} z ( i ) θ \theta θ
∂ ∂ θ [ ( z ( i ) + 34 ) log θ + 38 log ( 1 − θ ) ] = z ( i ) + 34 θ − 38 1 − θ = 0 \frac{\partial}{\partial\theta}\Big[(z^{(i)}+34)\log\theta + 38\log(1-\theta)\Big] = \frac{z^{(i)}+34}{\theta} - \frac{38}{1-\theta} = 0 ∂ θ ∂ [ ( z ( i ) + 34 ) log θ + 38 log ( 1 − θ ) ] = θ z ( i ) + 34 − 1 − θ 38 = 0
解得迭代公式
θ ( i + 1 ) = z ( i ) + 34 z ( i ) + 34 + 38 = z ( i ) + 34 z ( i ) + 72 \theta^{(i+1)} = \frac{z^{(i)} + 34}{z^{(i)} + 34 + 38} = \frac{z^{(i)} + 34}{z^{(i)} + 72} θ ( i + 1 ) = z ( i ) + 34 + 38 z ( i ) + 34 = z ( i ) + 72 z ( i ) + 34
反复用这两步迭代,θ \theta θ
假设观测数据 y 1 , y 2 , … , y N y_1, y_2, \dots, y_N y 1 , y 2 , … , y N
P ( y ∣ θ ) = ∑ k = 1 K α k ϕ ( y ∣ μ k , σ k 2 ) P(y\mid\theta) = \sum_{k=1}^{K} \alpha_k\, \phi(y\mid\mu_k,\sigma_k^2) P ( y ∣ θ ) = ∑ k = 1 K α k ϕ ( y ∣ μ k , σ k 2 )
其中 α k \alpha_k α k k k k α k ≥ 0 \alpha_k \ge 0 α k ≥ 0 ∑ k α k = 1 \sum_k \alpha_k = 1 ∑ k α k = 1 ϕ ( y ∣ μ k , σ k 2 ) \phi(y\mid\mu_k,\sigma_k^2) ϕ ( y ∣ μ k , σ k 2 ) k k k θ = { α k , μ k , σ k 2 } \theta = \{\alpha_k, \mu_k, \sigma_k^2\} θ = { α k , μ k , σ k 2 }
定义隐变量 γ j k \gamma_{jk} γ j k
γ j k = { 1 , 第 j 个观测来自第 k 个分模型 0 , 否则 \gamma_{jk} = \begin{cases} 1, & \text{第 } j \text{ 个观测来自第 } k \text{ 个分模型} \\ 0, & \text{否则} \end{cases} γ j k = { 1 , 0 , 第 j 个观测来自第 k 个分模型 否则
有了观测数据 y j y_j y j γ j k \gamma_{jk} γ j k 完全数据 就是 ( y j , γ j 1 , … , γ j K ) (y_j, \gamma_{j1}, \dots, \gamma_{jK}) ( y j , γ j 1 , … , γ j K )
P ( y , γ ∣ θ ) = ∏ k = 1 K α k n k ∏ j = 1 N [ ϕ ( y j ∣ μ k , σ k 2 ) ] γ j k , n k = ∑ j = 1 N γ j k P(y,\gamma\mid\theta) = \prod_{k=1}^{K} \alpha_k^{\,n_k} \prod_{j=1}^{N} \Big[\phi(y_j\mid\mu_k,\sigma_k^2)\Big]^{\gamma_{jk}}, \qquad n_k = \sum_{j=1}^{N}\gamma_{jk} P ( y , γ ∣ θ ) = ∏ k = 1 K α k n k ∏ j = 1 N [ ϕ ( y j ∣ μ k , σ k 2 ) ] γ j k , n k = ∑ j = 1 N γ j k
对应的完全数据对数似然 为
log P ( y , γ ∣ θ ) = ∑ k = 1 K { n k log α k + ∑ j = 1 N γ j k [ log 1 2 π − log σ k − 1 2 σ k 2 ( y j − μ k ) 2 ] } \log P(y,\gamma\mid\theta) = \sum_{k=1}^{K}\Big\{ n_k\log\alpha_k + \sum_{j=1}^{N}\gamma_{jk}\Big[\log\tfrac{1}{\sqrt{2\pi}} - \log\sigma_k - \tfrac{1}{2\sigma_k^2}(y_j-\mu_k)^2\Big]\Big\} log P ( y , γ ∣ θ ) = ∑ k = 1 K { n k log α k + ∑ j = 1 N γ j k [ log 2 π 1 − log σ k − 2 σ k 2 1 ( y j − μ k ) 2 ] }
E 步:确定 Q 函数 。把 γ j k \gamma_{jk} γ j k γ ^ j k \hat\gamma_{jk} γ ^ j k j j j k k k
γ ^ j k = α k ϕ ( y j ∣ μ k , σ k 2 ) ∑ l = 1 K α l ϕ ( y j ∣ μ l , σ l 2 ) \hat\gamma_{jk} = \frac{\alpha_k\,\phi(y_j\mid\mu_k,\sigma_k^2)}{\sum_{l=1}^{K}\alpha_l\,\phi(y_j\mid\mu_l,\sigma_l^2)} γ ^ j k = ∑ l = 1 K α l ϕ ( y j ∣ μ l , σ l 2 ) α k ϕ ( y j ∣ μ k , σ k 2 )
把 γ ^ j k \hat\gamma_{jk} γ ^ j k n k = ∑ j γ ^ j k n_k=\sum_j\hat\gamma_{jk} n k = ∑ j γ ^ j k
M 步 :对各参数求偏导并令其为 0,得到更新公式:
μ k = ∑ j = 1 N γ ^ j k y j ∑ j = 1 N γ ^ j k , σ k 2 = ∑ j = 1 N γ ^ j k ( y j − μ k ) 2 ∑ j = 1 N γ ^ j k , α k = ∑ j = 1 N γ ^ j k N = n k N \mu_k = \frac{\sum_{j=1}^{N}\hat\gamma_{jk}\,y_j}{\sum_{j=1}^{N}\hat\gamma_{jk}}, \qquad \sigma_k^2 = \frac{\sum_{j=1}^{N}\hat\gamma_{jk}(y_j-\mu_k)^2}{\sum_{j=1}^{N}\hat\gamma_{jk}}, \qquad \alpha_k = \frac{\sum_{j=1}^{N}\hat\gamma_{jk}}{N} = \frac{n_k}{N} μ k = ∑ j = 1 N γ ^ j k ∑ j = 1 N γ ^ j k y j , σ k 2 = ∑ j = 1 N γ ^ j k ∑ j = 1 N γ ^ j k ( y j − μ k ) 2 , α k = N ∑ j = 1 N γ ^ j k = N n k
E 步、M 步交替迭代到收敛,就得到 GMM 的参数估计。
在高斯混合模型下,当数据量 y 1 , … , y N y_1,\dots,y_N y 1 , … , y N
在 E 步 把数据分散到各个节点,每个节点在本地算出各自数据的条件概率(响应度) γ ^ j k \hat\gamma_{jk} γ ^ j k 期望充分统计量 (如 ∑ γ ^ j k \sum \hat\gamma_{jk} ∑ γ ^ j k ∑ γ ^ j k y j \sum \hat\gamma_{jk}y_j ∑ γ ^ j k y j ∑ γ ^ j k y j 2 \sum \hat\gamma_{jk}y_j^2 ∑ γ ^ j k y j 2
到 M 步 再把各节点的这些充分统计量汇总相加,统一更新全局参数 α k , μ k , σ k 2 \alpha_k, \mu_k, \sigma_k^2 α k , μ k , σ k 2
由于充分统计量是可加的,这种"各节点局部求和 + 全局汇总"的方式与单机 EM 结果完全一致,却能利用 Spark 这类框架并行处理海量数据。