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

Question

I am new to R as well as Bayesian Statistics. I am going through the problem set in Chapter#12 of Students Guide to Bayesian Statistics (this link has problem as well as answer plot).

In it Problem 12.4.3, author has provided a graph as the error vs number of samples.

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?

I am getting a straight line with no difference when sampling size is increased.

What am I missing? My R code:

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")
}

score 0 · Accepted Answer · answered Aug 09 '20 at 15:37

First. One issue is the way you compute the error, i.e. mean(replicate(1000, mean(rbinom(n=s, size = 1, prob=1/2 + ep)))) more or less equal to the probability of success, i.e. 1/2 + ep. Instead use the standard deviation or error.

Second. Besides this you are not really simulating a markovian coin, i.e. in your code the probabiliy of throwing "head" does not depend on the result of the last throw, e.g. rbinom(n=s, size = 1, prob=1/2 + ep) gives you s independent tosses of a coin with a probability of success of 1/2 + ep.

My code is probably not the most efficent solution to this kind of problem but it works. Try this:

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")

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

1 Answers1