MC, MCMC, Gibbs采样 原理&实现(in R)

Posted

tags:

篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了MC, MCMC, Gibbs采样 原理&实现(in R)相关的知识,希望对你有一定的参考价值。

本文用讲一下指定分布的随机抽样方法:MC(Monte Carlo), MC(Markov Chain), MCMC(Markov Chain Monte Carlo)的基本原理,并用R语言实现了几个例子:

1. Markov Chain (马尔科夫链)

2. Random Walk(随机游走)

3. MCMC具体方法:

     3.1 M-H法

     3.2 Gibbs采样 

 

PS:本篇blog为ese机器学习短期班参考资料(20140516课程),课上讲详述。

 

 

下面三节分别就前面几点简要介绍基本概念,并附上代码。这里的概念我会用最最naive的话去概括,详细内容就看我最下方推荐的链接吧(*^__^*) 

 

0. MC(Monte Carlo) 

     生成指定分布的随机数的抽样。

 

1. Markov Chain (马尔科夫链)

     假设 f(t) 是一个时间序列,Markov Chain是假设f(t+1)只与f(t)有关的随机过程。

     Implement in R:

 

 

[python] view plain copy
 
 技术分享技术分享
  1. #author: rachel @ ZJU  
  2. #email: [email protected]  
  3.   
  4. N = 10000  
  5. signal = vector(length = N)  
  6. signal[1] = 0  
  7. for (i in 2:N)  
  8. {  
  9.     # random select one offset (from [-1,1]) to signal[i-1]  
  10.     signal[i] = signal[i-1] + sample(c(-1,1),1)   
  11. }  
  12.   
  13. plot( signal,type = ‘l‘,col = ‘red‘)  



 

 

 

技术分享



2. Random Walk(随机游走)

如布朗运动,只是上面Markov Chain的二维拓展版:

Implement in R:

 

 

[python] view plain copy
 
 技术分享技术分享
  1. #author: rachel @ ZJU  
  2. #email: [email protected]  
  3.   
  4. N = 100  
  5. x = vector(length = N)  
  6. y = vector(length = N)  
  7. x[1] = 0  
  8. y[1] = 0  
  9. for (i in 2:N)  
  10. {  
  11.     x[i] = x[i-1] + rnorm(1)  
  12.     y[i] = y[i-1] + rnorm(1)  
  13. }  
  14.   
  15.   
  16. plot(x,y,type = ‘l‘, col=‘red‘)  



 

 

 

 

 

技术分享

 

 

 

3. MCMC具体方法:

    

MCMC方法最早由Metropolis(1954)给出,后来Metropolis的算法由Hastings改进,合称为M-H算法。M-H算法是MCMC的基础方法。由M-H算法演化出了许多新的抽样方法,包括目前在MCMC中最常用的Gibbs抽样也可以看做M-H算法的一个特例[2]。

 

概括起来,MCMC基于这样的理论,在满足【平衡方程】(detailed balance equation)条件下,MCMC可以通过很长的状态转移到达稳态。

【平衡方程】:
pi(x) * P(y|x) = pi(y) * P(x|y)
 
其中pi指分布,P指概率。这个平衡方程也就是表示条件概率(转化概率)与分布乘积的均衡.
 

 3.1 M-H法

1. 构造目标分布,初始化x0

2. 在第n步,从q(y|x_n) 生成新状态y

3. 以一定概率((pi(y) * P(x_n|y)) / (pi(x) * P(y|x_n)))接受y <PS: 看看上面的平衡方程,这个概率表示什么呢?参考这里[1]>

 

implementation in R:

 

 

[python] view plain copy
 
 技术分享技术分享
  1. #author: rachel @ ZJU  
  2. #email: [email protected]  
  3.   
  4. N = 10000  
  5. x = vector(length = N)  
  6. x[1] = 0  
  7.   
  8. # uniform variable: u  
  9. u = runif(N)  
  10. m_sd = 5  
  11. freedom = 5  
  12.   
  13. for (i in 2:N)  
  14. {  
  15.     y = rnorm(1,mean = x[i-1],sd = m_sd)  
  16.     print(y)  
  17.     #y = rt(1,df = freedom)  
  18.       
  19.     p_accept = dnorm(x[i-1],mean = y,sd = abs(2*y+1)) / dnorm(y, mean = x[i-1],sd = abs(2*x[i-1]+1))  
  20.     #print (p_accept)  
  21.       
  22.       
  23.     if ((u[i] <= p_accept))  
  24.     {  
  25.         x[i] = y  
  26.         print("accept")  
  27.     }  
  28.     else  
  29.     {  
  30.         x[i] = x[i-1]  
  31.         print("reject")  
  32.     }  
  33. }  
  34.   
  35. plot(x,type = ‘l‘)  
  36. dev.new()  
  37. hist(x)  



 

 


技术分享

 

 

技术分享

 

 

  3.2 Gibbs采样 

 
第n次,Draw 技术分享 from 技术分享,迭代采样结果接近真实p(\theta_1, \theta_2, ...)
 
也就是每一次都是固定其他参数,对一个参数进行采样。比如对于二元正态分布,其两个分量的一元条件分布仍满足正态分布:

技术分享
 
技术分享

 

那么在Gibbs采样中对其迭代采样的过程,实现如下:

 

 

[html] view plain copy
 
 技术分享技术分享
  1. #author: rachel @ ZJU  
  2. #email: [email protected]  
  3. #define Gauss Posterior Distribution  
  4.   
  5. p_ygivenx <- function(x,m1,m2,s1,s2)  
  6. {  
  7.     return (rnorm(1,m2+rho*s2/s1*(x-m1),sqrt(1-rho^2)*s2 ))  
  8. }  
  9.   
  10. p_xgiveny <- function(y,m1,m2,s1,s2)  
  11. {  
  12.     return  (rnorm(1,m1+rho*s1/s2*(y-m2),sqrt(1-rho^2)*s1 ))  
  13. }  
  14.   
  15.   
  16. N = 5000  
  17. K = 20 #iteration in each sampling  
  18. x_res = vector(length = N)  
  19. y_res = vector(length = N)  
  20. m1 = 10; m2 = -5; s1 = 5; s2 = 2  
  21. rho = 0.5  
  22. y = m2  
  23.   
  24. for (i in 1:N)  
  25. {  
  26.     for(i in 1:K)  
  27.     {  
  28.         x = p_xgiveny(y, m1,m2,s1,s2)  
  29.         y = p_ygivenx(x, m1,m2,s1,s2)  
  30.         # print(x)  
  31.         x_res[i] = x;  
  32.         y_res[i] = y;  
  33.     }  
  34. }  
  35.   
  36. hist(x_res,freq = 1)  
  37. dev.new()  
  38. plot(x_res,y_res)  
  39. library(MASS)  
  40. valid_range = seq(from = N/2, to = N, by = 1)  
  41. MVN.kdensity <- kde2d(x_res[valid_range], y_res[valid_range], h = 10) #估计核密度  
  42. plot(x_res[valid_range], y_res[valid_range], col = "blue", xlab = "x", ylab = "y")   
  43. contour(MVN.kdensity, add = TRUE)#二元正态分布等高线图  
  44.   
  45. #real distribution  
  46. real = mvrnorm(N,c(m1,m2),diag(c(s1,s2)))  
  47. # dev.new()  
  48. # plot(real[1:N,1],real[1:N,2])  



 

 

 

 


x分布图:

 

 

技术分享

 

 

(x,y)分布图:

技术分享

 

 

 

 

 

 

Reference:

1. http://www2.isye.gatech.edu/~brani/isyebayes/bank/handout10.pdf

2. http://site.douban.com/182577/widget/notes/10567181/note/292072927/

3. book:     http://statweb.stanford.edu/~owen/mc/

4. Classic: http://cis.temple.edu/~latecki/Courses/RobotFall07/PapersFall07/andrieu03introduction.pdf

from: http://blog.csdn.net/abcjennifer/article/details/25908495

以上是关于MC, MCMC, Gibbs采样 原理&实现(in R)的主要内容,如果未能解决你的问题,请参考以下文章

文本主题模型之LDA LDA求解之Gibbs采样算法

使用R语言进行Metroplis-in-Gibbs采样和MCMC运行分析

MCMC&Gibbs sampling

MCMCMCMC采样和M-H采样

MCMC蒙特卡罗方法

MCMC马尔科夫链