第19章 马尔可夫链蒙特卡洛法 - 二维Gibbs采样实例python实现 - 实现p_ygivenx, p_xgiveny函数

import numpy as np
import random

def p_ygivenx(x, mu1, mu2, sigma1, sigma2, rho):
    return np.random.normal(mu2 + rho * sigma2 / sigma1 * (x - mu1), np.sqrt(1 - rho**2) * sigma2)

def p_xgiveny(y, mu1, mu2, sigma1, sigma2, rho):
    return np.random.normal(mu1 + rho * sigma1 / sigma2 * (y - mu2), np.sqrt(1 - rho**2) * sigma1)

def gibbs_sampling(n_samples=1000):
    mu1, mu2 = 0, 0
    sigma1, sigma2 = 1, 1
    rho = 0.5
    
    samples = []
    x, y = 0, 0
    for _ in range(n_samples):
        x = p_xgiveny(y, mu1, mu2, sigma1, sigma2, rho)
        y = p_ygivenx(x, mu1, mu2, sigma1, sigma2, rho)
        samples.append((x, y))
    return np.array(samples)

samples = gibbs_sampling(100)
print("采样完成，均值:", np.mean(samples, axis=0))