计算机什么课学mcmc,MCMC例子

Love The Way You Lie 2022-10-07 13:47 114阅读 0赞

''''

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