本文用講一下指定分佈的隨機抽樣方法：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: [email protected]

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: [email protected]

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可以通過很長的狀態轉移到達穩態。

 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: [email protected]

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: [email protected]
#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,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, 後續內容繼續更新哦~