参考文章:
一文详解LDA主题模型
小白都能看懂的蒙特卡洛方法以及python实现
机器之心——马尔可夫模型
马尔可夫链蒙特卡罗算法(MCMC)
MCMC(一)蒙特卡罗方法
MCMC(二)马尔科夫链
MCMC(三)MCMC采样和M-H采样
通俗理解LDA模型
机器学习---算法基础(十二)LDA主题模型
数学基础
Gamma分布
Beta分布与Dirichlet分布
随机样本采样
蒙特卡洛方法
基本思想
蒙特卡罗方法也称统计模拟方法,是1940年代中期由于科学技术的发展和电子计算机的发明,而提出的一种以概率统计理论为指导的数值计算方法。是指使用随机数(或更常见的伪随机数)来解决很多计算问题的方法。
蒙特卡洛思想可以粗略的分成两类。
- 求解问题的内在随机性。借助计算机模拟这种随机过程
- 求解问题可以转化为某种随机分布的特征数,比如随机事件出现的概率,或者随机变量的期望值。通过随机抽样的方法,以随机事件出现的频率估计其概率,或者以抽样的数字特征估算随机变量的数字特征,并将其作为问题的解。这种方法多用于求解复杂的多维积分问题。
蒙特卡洛求积分
如果需要求一个分布在某个区间内的可能性,一般我们会采用如下的方式进行计算:
I = ∫ a b f ( x ) d x I=\int_a^b{f(x)dx} I=∫abf(x)dx
对于一个未知的分布,我们在不知道概率密度的情况系就很难计算积分。对于这种情况我们可以使用一种方式来计算某一种分布。
当我们在[a,b]之间随机取一点x时,它对应的函数值就是f(x)。接下来我们就可以用f(x) * (b - a)来粗略估计曲线下方的面积,也就是我们需要求的积分值,当然这种估计(或近似)是非常粗略的。
通过计算我们可以估计出,这个区间的概率分布。
S = 1 4 ( b − a ) ∑ 1 4 f ( x i ) S=\frac14(b-a)\sum_1^4{f(x_i)} S=41(b−a)1∑4f(xi)
当我们将这个概率上的分布取无限多个时,我们假设假设概率分布为 g ( x ) g(x) g(x),抽样的概率密度函数为f(x)且满足概率密度为 ∫ a b f ( x ) d x = 1 \int_a^bf(x)dx=1 ∫abf(x)dx=1(是在之前的抽样中其实是有一个假设的——g(x)在区间[a,b]上是连续分布的。但是很多情况下分布并不是连续的,所以取样采用的是概率分布的形式进行计算)。那么令 g ∗ ( x ) = g ( x ) f ( x ) g^\ast(x)=\frac{g(x)}{f(x)} g∗(x)=f(x)g(x),原有的积分可以写成如下形式:
I = ∫ a b g ∗ ( X ) f ( x ) d x I=\int_a^b g^\ast(X)f(x)dx I=∫abg∗(X)f(x)dx
由于采样时为离散的采样,所以计算公式可以修改为:
I = 1 N ∑ N g ∗ ( X ) f ( x ) I=\frac{1}{N}\sum_N g^\ast(X)f(x) I=N1N∑g∗(X)f(x)
一般取f(x)为均匀分布,当区间在[a,b]时所得到的积分为
I = b − a N ∑ N g ∗ ( X ) I=\frac{b-a}{N}\sum_N g^\ast(X) I=Nb−aN∑g∗(X)
拒绝接受采样(Acceptance-Rejection Sampling)
根据上一节我们可以知道,当我们了解样本的分布为常用分布时,我们可以通过蒙特卡洛方法来求解区间的概率值。
但是当我们无法了解样本的分布的时候我们就难以使用蒙克卡罗方法了。对于概率分布不是常见的分布,一个可行的办法是采用接受-拒绝采样来得到该分布的样本。既然 p ( x ) p(x) p(x) 太复杂在程序中没法直接采样,那么我设定一个程序可采样的分布 q ( x ) q(x) q(x)比如高斯分布,然后按照一定的方法拒绝某些样本,以达到接近 p ( x ) p(x) p(x) 分布的目的,其中 q ( x ) q(x) q(x)叫做 proposal distribution(拒绝接收采样)
具体采用过程如下,设定一个方便采样的常用概率分布函数 q ( x ) q(x) q(x),以及一个常量 k k k,使得 p ( x ) p(x) p(x) 总在 k q ( x ) k_q(x) kq(x) 的下方。下图所示。
首先,采样得到 q ( x ) q(x) q(x)的一个样本 z 0 z_0 z0,采样方法如第三节。然后,从均匀分布 ( 0 , k q ( z 0 ) ) (0,k_q(z_0)) (0,kq(z0))中采样得到一个值 u u u。如果 u u u落在了上图中的灰色区域,则拒绝这次抽样,否则接受这个样本 z 0 z_0 z0。重复以上过程得到n个接受的样本 z 0 , z 1 , . . . z n − 1 z_0,z_1,...z_{n−1} z0,z1,...zn−1,则最后的蒙特卡罗方法求解结果为:
1 n ∑ i = 0 n − 1 f ( z i ) p ( z i ) \frac1n\sum_{i=0}^{n-1}\frac{f(z_i)}{p(z_i)} n1i=0∑n−1p(zi)f(zi)
通过拒绝采样和蒙特卡洛方法我们依然不容易求出所需要的样本点。
- 对于一些二维分布p(x,y),有时候我们只能得到条件分布p(x|y)和p(y|x)和,却很难得到二维分布p(x,y)一般形式,这时我们无法用接受-拒绝采样得到其样本集。
- 对于一些高维的复杂非常见分布p(x1,x2,…,xn),我们要找到一个合适的q(x)和k非常困难。
通过上述可以看出,对于复杂概率模型,蒙特卡洛与拒绝采样方法难以应付,此时需要请出马尔科夫模型来对复杂概率模型进行样本的采样。
马尔科夫模型
马尔科夫性质
马尔可夫模型是一类随机过程,该过程有如下的特点:已知有几种状态,各种状态之间可以互相转化。 状态空间中经过从一个状态到另一个状态的转换的随机过程。该过程要求具备“无记忆”的性质:下一状态的概率分布只能由当前状态决定,在时间序列中它前面的事件均与之无关。 这种特定类型的“无记忆性”称作马尔可夫性质。马尔科夫链作为实际过程的统计模型具有许多应用。其公式如下:
模型在t+1时刻的状态条件概率仅仅依赖于t时刻。
马尔科夫模型
「马尔可夫模型」是指基于马尔可夫性质的模型,其假设一个给定过程的未来状态仅取决于当前状态。根据系统状态是否完全可被观测以及系统是自动的还是受控的,可以将常见的马尔可夫模型分成四种:马尔可夫链、隐马尔可夫模型(HMM)、马尔可夫决策过程(MDP)和部分可观测马尔可夫决策过程(POMDP),具体可见下表。
最简单的马尔可夫模型是马尔可夫链。它用一个随时间变化的随机变量来模拟系统的状态。在这种情况下,马尔可夫性质表明,这个变量的分布只取决于之前状态的分布。当马尔可夫链的状态只能部分观察到,这就是隐马尔可夫模型。隐马尔可夫模型常用的用途是语音识别,它是大多数现代自动语音识别系统的基础。
建立一个马尔科夫模型如下所示:
其状态之间的转移我们用转移概率矩阵P表示:
此时,我们给定一个初始状态,然后经过该状态转移矩阵的转换,最终会收敛到一个稳定的状态,具体如马尔科夫链定理所示
当马尔科夫链达到平稳分布的情况下,我们得到的性质为:从某个时刻开始,每个状态之间发生的概率不再变化,我们成这个最后不在变化的概率为马尔科夫链的平稳分布。
由于马尔科夫链能收敛到平稳分布, 于是有了一个想法:如果我们能构造一个转移矩阵为P的马氏链,使得该马氏链的平稳分布恰好是p(x), 那么我们从任何一个初始状态x0出发沿着马氏链转移, 得到一个转移序列 x0, x1, x2,⋯xn, xn+1⋯, 如果马氏链在第n步已经收敛了,于是我们就得到了 π(x) 的样本xn, xn+1⋯(也就是从第n步收敛时开始,之后的x都服从同一个平稳分布,我们可以将这个分布设定为我们的目标采样分布)。
从上面可以看出马尔科夫链的平稳分布收敛主要依赖于状态转移矩阵,所以关键是如何构建状态转移矩阵,使得最终的平稳分布是我们所要的分布。想做到这一点主要依赖于细致平稳定理
马尔可夫采样
如果我们得到了某个平稳分布所对应的马尔科夫链状态转移矩阵,我们就很容易采用出这个平稳分布的样本集。
假设我们任意初始的概率分布是 π 0 ( x ) π_0(x) π0(x), 经过第一轮马尔科夫链状态转移后的概率分布是 π 1 ( x ) π_1(x) π1(x),。。。第i轮的概率分布是 π i ( x ) π_i(x) πi(x)。假设从第i轮开始再经过n轮后马尔科夫链收敛到我们的平稳分布 π ( x ) π(x) π(x),即:
π n ( x ) = π n + 1 ( x ) = π n + 2 ( x ) = . . . = π ( x ) π_n(x)=π_{n+1}(x)=π_{n+2}(x)=...=π(x) πn(x)=πn+1(x)=πn+2(x)=...=π(x)
对于每个分布 π i ( x ) π_i(x) πi(x),我们有:
π i ( x ) = π i − 1 ( x ) P = π i − 2 ( x ) P 2 = π 0 ( x ) P i π_i(x)=π_{i−1}(x)P=π_{i−2}(x)P^2=π_0(x)P^i πi(x)=πi−1(x)P=πi−2(x)P2=π0(x)P