MC, MCMC, Gibbs Sampling principle & implementation (in R)

This article uses a random sampling method that specifies the distribution: MC (Monte Carlo), MC (Markov Chain), MCMC (Markov Chain Monte Carlo), and uses the R language to implement several examples:

1. Markov Chain (Markov chain)

2. Random Walk (stochastic walk)

3. MCMC Specific methods:

3.1 M-h method

3.2 Gibbs Sampling

PS: This blog is an ESE machine learning short class reference (20140516 course), detailed in class.

The following sections briefly describe the basic concepts and attach the code to the preceding points. Here the concept I will use the most naive words to summarize, the details on the most recommended link to my bottom (*^__^*)

0. MC (Monte Carlo)

Generates a sampling of random numbers for a specified distribution.

1. Markov Chain (Markov chain)

Assuming that f (t) is a time series, Markov chain is a stochastic process that assumes that F (t+1) is only related to f (t).

Implement in R:

[Python]View PlainCopy 

2. #author: Rachel @ zju
4. #email: [Email protected]
8. N = 10000
10. signal = vector (length = N)
12. signal[1] = 0
14. For (i in 2:n)
16. {
18. # Random Select one offset (from [ -1,1]) to Signal[i-1]
20. Signal[i] = signal[i-1] + sample (C (-1,1),1)
22. }
26. Plot (Signal,type = ' l ', col = ' red ')

2. Random Walk (stochastic walk)

such as Brownian motion, just above Markov chain's two-dimensional expansion version:

Implement in R:

[Python]View PlainCopy 

2. #author: Rachel @ zju
4. #email: [Email protected]
8. N =
10. x = vector (length = N)
12. y = vector (length = N)
14. x[1] = 0
16. y[1] = 0
18. For (i in 2:n)
20. {
22. X[i] = x[i-1] + rnorm (1)
24. Y[i] = y[i-1] + rnorm (1)
26. }
32. Plot (X,y,type = ' l ', col=' red ')

3. MCMC Specific methods:

The MCMC method was first given by Metropolis (1954), and later the metropolis algorithm was improved by Hastings, collectively known as the M-H algorithm. The M-H algorithm is the basic method of MCMC. Many new sampling methods have been developed by the M-H algorithm, including the most commonly used Gibbs sampling in MCMC, which can be seen as a special case of the M-H algorithm [2].

To sum up, MCMC based on this theory, in satisfying the "equilibrium equation" (Detailed balance equation) conditions, the MCMC can be transferred to the steady state through a very long transition.

"Equilibrium equation": Pi (x) * p (y|x) = Pi (y) * p (x|y) where Pi refers to distribution, P-refers to probability. This equation of equilibrium is the equilibrium that represents the probability of the condition (conversion probability) and the product of the distribution.

3.1 M-h method

1. Construct the target distribution, initialize the x0

2. In nth step, generate a new state from Q (y|x_n) y

3. Accept Y <ps with a certain probability ((Pi (y) * p (X_n|y))/(PI (x) * p (y|x_n)): Look at the equilibrium equation above, what does this probability mean? Refer here and [1]>

Implementation in R:

[Python]View PlainCopy 

2. #author: Rachel @ zju
4. #email: [Email protected]
8. N = 10000
10. x = vector (length = N)
12. x[1] = 0
16. # Uniform Variable:u
18. U = runif (N)
20. M_SD = 5
22. Freedom = 5
26. For (i in 2:n)
28. {
30. y = rnorm (1,mean = x[i-1],sd = m_sd)
32. print (y)
34. #y = RT (1,DF = Freedom)
38. 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 ))  
40. #print (p_accept)
46. if ((U[i] <= p_accept))
48. {
50. X[i] = y
52. Print ("accept")
54. }
56. Else
58. {
60. X[i] = x[i-1]
62. print ("reject")
64. }
66. }
70. Plot (X,type = ' l ')
72. Dev.new ()
74. hist (x)

3.2 Gibbs Sampling

Nth time, Draw from, iterative sampling results close to real P (\theta_1, \theta_2, ...) that is, each time the other parameters are fixed, a parameter is sampled. For example, for a two-dollar normal distribution, the one-dimensional conditional distribution of its two components still satisfies the normal distribution:

Then the process of iterative sampling in Gibbs sampling is implemented as follows:

[HTML]View PlainCopy  

2. #author: Rachel @ zju
4. #email: [Email protected]
6. #define GAUSS Posterior distribution
10. P_ygivenx <-function (X,M1,M2,S1,S2)
12. {
14. Return (Rnorm (1,m2+rho*s2/s1* (X-M1), sqrt (1-rho^2) *s2))
16. }
20. P_xgiveny <-function (Y,M1,M2,S1,S2)
22. {
24. Return (Rnorm (1,m1+rho*s1/s2* (y-m2), sqrt (1-rho^2) *s1))
26. }
32. N =
34. K = #iteration in each sampling
36. X_res = vector (length = N)
38. Y_res = vector (length = N)
40. M1 = 10; m2 =-5; S1 = 5; s2 = 2
42. Rho = 0.5
44. y = m2
48. For (i in 1:n)
50. {
52. For (i in 1:k)
54. {
56. x = P_xgiveny (y, m1,m2,s1,s2)
58. y = p_ygivenx (x, m1,m2,s1,s2)
60. # Print (x)
62. X_res[i] = x;
64. Y_res[i] = y;
66. }
68. }
72. hist (x_res,freq = 1)
74. Dev.new ()
76. Plot (X_res,y_res)
78. Library (MASS)
80. Valid_range = seq (from= N/2, to = N, by = 1)
82. Mvn.kdensity <-Kde2d (X_res[valid_range], Y_res[valid_range], h = ten) #估计核密度
84. Plot (X_res[valid_range], Y_res[valid_range], col = "Blue", Xlab = "x", Ylab = "y")
86. Contour (mvn.kdensity, add = TRUE) #二元正态分布等高线图
90. #real distribution
92. # real = Mvrnorm (n,c (m1,m2), Diag (C (S1,S2)))
94. # dev.new ()
96. # Plot (real[1:n,1],real[1:n,2])

X Distribution Chart:

(x, y) distribution map:

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 Sampling principle & implementation (in R)