R语言用狄里克莱分布作为建议值的例子

以下是一个使用狄利克雷分布作为建议分布进行MCMC采样的R语言代码示例：

library(mvtnorm)  # 加载mvtnorm包，用于生成多元正态分布随机数

# 真实参数值
theta <- c(0.3, 0.4, 0.3)

# 生成数据集
n <- 1000
x <- rmultinom(n, size = 1, prob = theta)

# 定义目标函数
loglik <- function(theta, x) {
  sum(dmultinom(x, size = 1, prob = theta, log = TRUE))
}

# 定义狄利克雷分布的参数
alpha <- c(1, 1, 1)

# 定义MCMC采样函数
mcmc <- function(x, alpha, niter, burnin) {
  # 初始化参数值
  theta <- rep(0, length(alpha))
  theta[1] <- runif(1)
  theta[2:length(alpha)] <- rdirichlet(1, alpha[-1])
  
  # 初始化接受率
  accept <- 0
  
  # 初始化保存参数值的矩阵
  theta.mat <- matrix(0, nrow = niter, ncol = length(alpha))
  
  # 进行MCMC采样
  for (i in 1:niter) {
    # 更新theta[1]
    prop.theta.1 <- rnorm(1, theta[1], 0.1)
    prop.theta <- c(prop.theta.1, theta[-1])
    prop.loglik <- loglik(prop.theta, x)
    curr.loglik <- loglik(theta, x)
    alpha.ratio <- sum(alpha)
    prop.ratio <- sum(prop.theta)
    log.ratio <- alpha.ratio * log(prop.ratio) - alpha.ratio * log(sum(theta)) +
      sum((prop.theta - theta) * (alpha - 1) / theta) + prop.loglik - curr.loglik
    if (log(runif(1)) < log.ratio) {
      theta[1] <- prop.theta.1
      accept <- accept + 1
    }
    
    # 更新theta[2:length(alpha)]
    for (j in 2:length(alpha)) {
      prop.theta <- theta
      prop.theta[j] <- rnorm(1, theta[j], 0.1)
      prop.theta <- prop.theta / sum(prop.theta)
      prop.loglik <- loglik(prop.theta, x)
      curr.loglik <- loglik(theta, x)
      log.ratio <- (alpha[j] - 1) * log(prop.theta[j] / theta[j]) + prop.loglik - curr.loglik
      if (log(runif(1)) < log.ratio) {
        theta[j] <- prop.theta[j]
        accept <- accept + 1
      }
    }
    
    # 保存参数值
    theta.mat[i, ] <- theta
  }
  
  # 计算接受率
  accept.rate <- accept / (niter * (length(alpha) - 1))
  
  # 返回结果
  result <- list(theta.mat = theta.mat[(burnin + 1):niter, ], accept.rate = accept.rate)
  return(result)
}

# 进行MCMC采样
set.seed(123)
result <- mcmc(x, alpha, niter = 5000, burnin = 1000)

# 输出结果
summary(result$theta.mat)
plot(result$theta.mat[, 1], type = "l", xlab = "Iteration", ylab = expression(theta[1]))
plot(result$theta.mat[, 2], type = "l", xlab = "Iteration", ylab = expression(theta[2]))
plot(result$theta.mat[, 3], type = "l", xlab = "Iteration", ylab = expression(theta[3]))

这个例子中，我们使用了狄利克雷分布作为建议分布，并且使用了自适应MCMC算法来进行参数估计。在MCMC采样过程中，我们使用了随机游走元算法来更新参数值，以期望能够在参数空间中均匀地采样。最终，我们得到了参数的后验分布估计结果，并进行了可视化展示。