需要帮助 w/R 代码绘图错误与否。马尔可夫链中的样本

Question

我是 R 和贝叶斯统计的新手。我正在研究 Students Guide to Bayesian Statistics 第 12 章中设置的问题（这个 link 有问题也有答案）。

在问题 12.4.3 中，作者提供了一个误差与样本数的关系图。

Consider a type of coin for which the result of the next throw (heads or tails) can depend on the result of the current throw. In particular, if a heads is thrown then the probability of obtaining a heads on the next throw is ; if instead a tails is thrown then the probability of obtaining a tails on the next throw is . To start, we assume 0 ≤ ∈ ≤. The random variable X takes the value 0 if the coin lands tails up or 1 if it lands heads up on a given throw.

Problem 12.4.3 As ∈ increases, how does the error in estimating the mean change, and why?

当样本量增加时，我得到一条没有差异的直线。

我错过了什么？我的 R 代码：

epsilon <- seq(from = 0, to = 0.5, length.out = 10 )
first_throw <- rbinom(n=1, size=1, prob = 1/2)
cat("\nFirst Throw: ",first_throw)
last_throw <- first_throw


for ( s in c(10,20,100)){

  for (ep in epsilon) {
    j <- 1
    curr_err <- 0
      if(last_throw == 1){
        last_throw <- rbinom(n=1, size = 1,  prob=1/2 + ep)
        curr_err <- abs(mean(replicate(1000, mean(rbinom(n=s, size = 1,  prob=1/2 + ep)))) - 0.5)
    }
    else{
      last_throw <- rbinom(n=1, size = 1,  prob=1/2 - ep)
      curr_err <- abs(mean(replicate(1000, mean(rbinom(n=s, size = 1,  prob=1/2 - ep)))) - 0.5)
    }
    
    lerrors [j] <- curr_err
    j <- j + 1
    
  }
  cat("\n epsilon: ", epsilon)
  cat("\n lerrors: ", lerrors)
  plot(epsilon,lerrors, col="blue")
  lines(epsilon, lerrors, col="blue")
}

Answer 1

首先。一个问题是计算误差的方式，即 mean(replicate(1000, mean(rbinom(n=s, size = 1, prob=1/2 + ep)))) 或多或少等于成功概率，即 1/2 + ep。而是使用标准差或误差。

其次。除此之外，您并不是真正在模拟马尔可夫硬币，即在您的代码中，抛出“正面”的概率并不取决于最后一次抛出的结果，例如rbinom(n=s, size = 1, prob=1/2 + ep) 给你 s 次独立抛硬币，成功概率为 1/2 + ep。

我的代码可能不是解决此类问题的最有效方法，但它确实有效。试试这个：

set.seed(42)

markov_coin <- function(n, prob_head = .5, eps = .1) {
  res <- vector("integer", n)
  for (s in seq(n)) {
    res[s] <- if (s == 1) {
      rbinom(1, 1, prob_head)
    } else {
      rbinom(1, 1, prob_head + ifelse(res[s-1] == 1, eps, -eps))
    }
  }
  res
}

epsilon <- seq(0, .5, length.out = 10)
svec <- c(10, 20, 100)

res <- matrix(nrow = 10, ncol = 3)
for (s in seq_along(svec)) {
  for (eps in seq_along(epsilon)) {
    res[eps, s] <- sd(replicate(1000, mean(markov_coin(svec[[s]], .5, epsilon[eps])))) 
  }
}

plot(epsilon, res[,1], "b", col="blue", xlab = "epsilon", ylab = "error", ylim = c(0, .5))
lines(epsilon, res[,2], "b", col="orange")
lines(epsilon, res[,3], "b", col="green")

需要帮助 w/R 代码绘图错误与否。马尔可夫链中的样本

Need help w/ R code plotting error vs no. of samples in Markovian chain

r

bayesian