''''
假设目标分布是: 一维的正太分布, 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()