R语言 Gibbs sampler

Gibbs采样是一种马尔科夫链蒙特卡罗（MCMC）方法，常用于从联合分布中采样，特别是在贝叶斯统计中。在R语言中，可以使用Gibbs采样生成随机样本。假设我们想从二元正态分布中采样（X1，X2），其中X1和X2的条件分布为：

X1|X2 ~ N(μ1+ρσ1(X2-μ2), (1-ρ^2)σ1^2)
X2|X1 ~ N(μ2+ρσ2(X1-μ1), (1-ρ^2)σ2^2)

那么我们可以通过以下代码实现Gibbs采样：

gibbs <- function(N, mu1, mu2, sigma1, sigma2, rho) {
  x <- y <- numeric(N)
  x[1] <- y[1] <- 0
  for (i in 2:N) {
    x[i] <- rnorm(1, mean = mu1 + rho * sigma1 * (y[i-1] - mu2), sd = sqrt(1 - rho^2) * sigma1)
    y[i] <- rnorm(1, mean = mu2 + rho * sigma2 * (x[i] - mu1), sd = sqrt(1 - rho^2) * sigma2)
  }
  return(data.frame(x, y))
}

# 示例
set.seed(123)
samples <- gibbs(N = 1000, mu1 = 0, mu2 = 0, sigma1 = 1, sigma2 = 1, rho = 0.5)
plot(samples, type = "l", main = "Gibbs sampler for bivariate normal distribution")

这里我们使用了`rnorm()`函数生成正态分布随机数，其中mean参数为均值，sd参数为标准差。在每次迭代中，我们分别生成X1和X2的下一个样本值，并将其存储在向量x和y中。最后，我们返回一个数据框，其中包含了所有采样值，并将其可视化。