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

本文深入浅出地介绍了Monte Carlo方法、马尔科夫链(MarkovChain)及Markov Chain Monte Carlo(MCMC)方法的基本原理,并通过R语言实现了一系列实例,包括MarkovChain、随机游走(RandomWalk)和MCMC的具体应用。重点讲解了M-H算法和Gibbs采样,同时提供了实际代码实现,适合机器学习短期班的学习者参考。

摘要生成于 C知道 ,由 DeepSeek-R1 满血版支持, 前往体验 >

本文用讲一下指定分布的随机抽样方法: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:


#author: rachel @ ZJU
#email: zrqjennifer@gmail.com

N = 10000
signal = vector(length = N)
signal[1] = 0
for (i in 2:N)
{
	# random select one offset (from [-1,1]) to signal[i-1]
	signal[i] = signal[i-1] + sample(c(-1,1),1) 
}

plot( signal,type = 'l',col = 'red')






2. Random Walk(随机游走)

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

Implement in R:


#author: rachel @ ZJU
#email: zrqjennifer@gmail.com

N = 100
x = vector(length = N)
y = vector(length = N)
x[1] = 0
y[1] = 0
for (i in 2:N)
{
	x[i] = x[i-1] + rnorm(1)
	y[i] = y[i-1] + rnorm(1)
}


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:


#author: rachel @ ZJU
#email: zrqjennifer@gmail.com

N = 10000
x = vector(length = N)
x[1] = 0

# uniform variable: u
u = runif(N)
m_sd = 5
freedom = 5

for (i in 2:N)
{
	y = rnorm(1,mean = x[i-1],sd = m_sd)
	print(y)
	#y = rt(1,df = freedom)
	
	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))
	#print (p_accept)
	
	
	if ((u[i] <= p_accept))
	{
		x[i] = y
		print("accept")
	}
	else
	{
		x[i] = x[i-1]
		print("reject")
	}
}

plot(x,type = 'l')
dev.new()
hist(x)









  3.2 Gibbs采样 


第n次,Draw from ,迭代采样结果接近真实p(\theta_1, \theta_2, ...)

也就是每一次都是固定其他参数,对一个参数进行采样。比如对于二元正态分布,其两个分量的一元条件分布仍满足正态分布:



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


#author: rachel @ ZJU
#email: zrqjennifer@gmail.com
#define Gauss Posterior Distribution

p_ygivenx <- function(x,m1,m2,s1,s2)
{
	return (rnorm(1,m2+rho*s2/s1*(x-m1),sqrt(1-rho^2)*s2 ))
}

p_xgiveny <- function(y,m1,m2,s1,s2)
{
	return 	(rnorm(1,m1+rho*s1/s2*(y-m2),sqrt(1-rho^2)*s1 ))
}


N = 5000
K = 20 #iteration in each sampling
x_res = vector(length = N)
y_res = vector(length = N)
m1 = 10; m2 = -5; s1 = 5; s2 = 2
rho = 0.5
y = m2

for (i in 1:N)
{
	for(i in 1:K)
	{
		x = p_xgiveny(y, m1,m2,s1,s2)
		y = p_ygivenx(x, m1,m2,s1,s2)
		# print(x)
		x_res[i] = x;
		y_res[i] = y;
	}
}

hist(x_res,freq = 1)
dev.new()
plot(x_res,y_res)
library(MASS)
valid_range = seq(from = N/2, to = N, by = 1)
MVN.kdensity <- kde2d(x_res[valid_range], y_res[valid_range], h = 10) #估计核密度
plot(x_res[valid_range], y_res[valid_range], col = "blue", xlab = "x", ylab = "y") 
contour(MVN.kdensity, add = TRUE)#二元正态分布等高线图

#real distribution
# real = mvrnorm(N,c(m1,m2),diag(c(s1,s2)))
# dev.new()
# 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





欢迎参与讨论并关注本博客和微博Rachel____Zhang, 后续内容继续更新哦~





### Gibbs Sampling 算法原理实现 #### 1. 算法原理 Gibbs Sampling 是一种 Markov Chain Monte Carlo (MCMC) 方法,用于从复杂的联合分布中抽样。它通过在高维空间中逐步更新每个变量的值,从而生成样本序列。该算法的核心思想是基于条件分布进行采样,而不是直接从联合分布中采样[^1]。具体来说,假设有一个多维随机变量 \( \mathbf{X} = (X_1, X_2, \ldots, X_m) \),其联合分布为 \( P(X_1, X_2, \ldots, X_m) \)Gibbs Sampling 的步骤如下: - 初始化所有变量的值。 - 按照一定的顺序(例如从 \( X_1 \) 到 \( X_m \)),依次从条件分布 \( P(X_i | X_1, X_2, \ldots, X_{i-1}, X_{i+1}, \ldots, X_m) \)采样新的值。 - 重复上述过程,直到马尔可夫链达到平稳分布。 这种方法的优点在于,即使联合分布难以直接采样,只要条件分布容易计算,就可以有效地生成样本[^2]。 #### 2. 实现方法 以下是使用 Python 实现 Gibbs Sampling 的一个简单示例,假设目标分布为二维正态分布: ```python import numpy as np import matplotlib.pyplot as plt # 条件分布的定义 def p_x_given_y(y, mu_x=0, sigma_x=1, rho=0.8): mu = mu_x + rho * sigma_x * y sigma = np.sqrt(1 - rho**2) * sigma_x return np.random.normal(mu, sigma) def p_y_given_x(x, mu_y=0, sigma_y=1, rho=0.8): mu = mu_y + rho * sigma_y * x sigma = np.sqrt(1 - rho**2) * sigma_y return np.random.normal(mu, sigma) # Gibbs Sampling 过程 def gibbs_sampling(num_samples=1000, burn_in=500): samples = np.zeros((num_samples, 2)) x, y = 0, 0 # 初始化 for i in range(num_samples + burn_in): x = p_x_given_y(y) y = p_y_given_x(x) if i &gt;= burn_in: # 跳过 burn-in 阶段 samples[i - burn_in] = [x, y] return samples # 生成样本并可视化 samples = gibbs_sampling() plt.scatter(samples[:, 0], samples[:, 1], alpha=0.5, s=1) plt.title(&quot;Gibbs Sampling Results&quot;) plt.xlabel(&quot;X&quot;) plt.ylabel(&quot;Y&quot;) plt.show() ``` #### 3. 理论基础与优缺点 Gibbs Sampling 的理论基础建立在马尔可夫链的收敛性上。根据马尔可夫链的性质,经过足够多次迭代后,样本序列将趋于平稳分布。然而,Gibbs Sampling 也存在一些局限性: - 如果条件分布难以计算或采样,则该方法可能不适用。 - 在高维空间中,马尔可夫链可能存在混合缓慢的问题,导致收敛速度较慢[^3]。 #### 4. 应用场景 Gibbs Sampling 广泛应用于统计推断、贝叶斯学习和隐变量模型等领域。例如,在 LDA(Latent Dirichlet Allocation)模型中,Gibbs Sampling 用于估计话题分配 \( z \) 和话题-词分布 \( \phi \)[^1]。 ---
评论 28
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值