吉布斯采样python代码
时间: 2023-11-16 17:56:13 浏览: 43
以下是吉布斯采样的Python代码:
```python
import random
import math
import matplotlib.pyplot as plt
def p_ygivenx(x, m1, m2, s1, s2):
return (random.normalvariate(m2 + rho * s2 / s1 * (x - m1), math.sqrt(1 - rho ** 2) * s2))
def p_xgiveny(y, m1, m2, s1, s2):
return (random.normalvariate(m1 + rho * s1 / s2 * (y - m2), math.sqrt(1 - rho ** 2) * s1))
N = 5000
K = 20
x_res = []
y_res = []
m1 = 10
m2 = -5
s1 = 5
s2 = 2
rho = 0.5
y = m2
for i in range(N):
for j in range(K):
x = p_xgiveny(y, m1, m2, s1, s2)
y = p_ygivenx(x, m1, m2, s1, s2)
x_res.append(x)
y_res.append(y)
num_bins = 50
plt.hist(x_res, num_bins, normed=1, facecolor='green', alpha=0.5)
plt.hist(y_res, num_bins, normed=1, facecolor='red', alpha=0.5)
plt.title('Histogram')
plt.show()
```
这段代码实现了一个简单的吉布斯采样过程,用于估计分布的期望和边缘分布。其中,`p_ygivenx`和`p_xgiveny`分别表示给定$x$时$y$的条件概率和给定$y$时$x$的条件概率,$N$表示采样次数,$K$表示每次采样的迭代次数,$x_res$和$y_res$分别表示采样得到的$x$和$y$的结果。最后,代码使用`matplotlib`库绘制了$x$和$y$的直方图。