''''

假设目标分布是: 一维的正太分布, i.e., N(0, 1)

建议的转移矩阵Q(i,j)也是正太分布， j 服从 N(i, 2^2)

pi(i)Q(i,j)*alpha(i,j) = pi(j)Q(i,j)*alpha(i,j)

where alpha(i,j) = pi(j)Q(j,i) 接受率

alpha(i,j) 可能很小，导致接受率比较小，采样效率低下

修正: alpha(i,j) = min(1, alpha(i,j ) / alpha(j, i))

= min(1, pi(j)Q(j,i) / pi(i)Q(i,j))

= min(1, pi(j) / pi(i)) # 通常Q(i,j)是对称的, 并且此时分布pi都无需归一化

Gibbs采样：

考虑二维的情况, A = (x1, y1), B = (x1, y2), C = (x2, y2)

那么有： pi(A)pi(B|A) = p(x1)p(y1|x1) * p(y2|x1)

pi(B)pi(A|B) = p(x1)p(y2|x1) * p(y1|x1)

上面两个式子相等：

==> pi(A)*p(y2|x1) = pi(B)*p(y1|x1)， A和B具有相同的x坐标x1

同理，我们有

pi(B)*p(x2|y2) = pi(C)*p(x1|y2)， B和C具有相同的y坐标y2

因此，我们可以构造出接受率为1的状态转移矩阵，如上定义，沿着某个坐标轴进行采样，其余坐标保持不变

高维的情况类似'''

from scipy importstatsimportnumpy as npimportmatplotlib.pyplot as pltdefMCMC_01():deftarget_pdf(x):return stats.norm.pdf(x, loc=0, scale=1) #随机变量的概率密度函数

X_0= 2.0Sample_X=[X_0]for t in range(100000):

X_i= Sample_X[-1]

X_j= np.random.randn() * 2 + X_i #服从分布 N(i, 2)

alpha = min(1, target_pdf(x=X_j) / target_pdf(x=X_i))

p=np.random.rand()if p < alpha: #接受

Sample_X.append(X_j)else:

Sample_X.append(X_i)

show_cnt= 10000Sample_X= Sample_X[-show_cnt:]

Standard_X=np.random.randn(show_cnt)

plt.figure(figsize=(15, 9))

plt.subplot(1, 2, 1)

plt.hist(Sample_X, bins=50, normed=1)

plt.xlabel('x')

plt.xlabel('p')

plt.title('MCMC')

plt.subplot(1, 2, 2)

plt.hist(Standard_X, bins=50, normed=1)

plt.xlabel('x')

plt.title('Normal(0, 1)')

plt.show()

MCMC_01()defMCMC_02():'''目标分布是二维的正太分布：N (u, sigma^2),

u= (u1, u2) = (3, 5),

sigma = [[sigma_1^2, rho*sigma_1*sigma_2],

[rho*sigma_1*sigma_2, sigma_2^2]]

= [[1, 2],

[2, 4]]

那么对应的条件分布为：

P(x1|x2) = N(u1 + (x2 - u2) * rho * sigma_1 / sigma_2, (1-rho^2) * sigma_1^2)

P(x2|x1) = N(u2 + (x1 - u1) * rho * sigma_2 / sigma_1, (1-rho^2) * sigma_2^2)'''u1, u2= 3, 5sigma_1, sigma_2, rho= 1, 2, 0.5

defgibbs_sampling_x2_given_x1(x1):#return: p(x2|x1)

now_u = u2 + (x1 - u1) * rho * sigma_2 /sigma_1

now_sigma= np.sqrt(1 - rho ** 2) * sigma_2 ** 2now_sigma=np.sqrt(now_sigma)return np.random.randn() * now_sigma +now_udefgibbs_sampling_x1_given_x2(x2):#return: p(x1|x2)

now_u = u1 + (x2 - u2) * rho * sigma_1 /sigma_2

now_sigma= (1-rho ** 2) * sigma_1 ** 2now_sigma=np.sqrt(now_sigma)return np.random.randn() * now_sigma +now_u

X_0= [-2.0, -3.0]

Sample_X=[X_0]for t in range(100000):

x_2_pre= Sample_X[-1][-1]

x_1= gibbs_sampling_x1_given_x2(x2=x_2_pre) #根据上一个采样点的x2, 采样现在的x1

x_2 = gibbs_sampling_x2_given_x1(x1=x_1) #根据刚刚采样点的x1, 采样现在的x2

Sample_X.append([x_1, x_2])if t<10:print('x_1:{} x_2:{}'.format(x_1, x_2))

show_cnt= 10000Sample_X= np.array(Sample_X[-show_cnt:])#Standard_X = np.random.randn(show_cnt)

print('Sample_X:', Sample_X.shape)

plt.figure(figsize=(10, 9))

plt.subplot(1, 1, 1)

plt.plot(Sample_X[:, 0], Sample_X[:,1], 'r.')

plt.xlabel('$x_1$')

plt.ylabel('$x_2$')

plt.title('Gibbs')

plt.show()

MCMC_02()