Chapter 7, Exercise 15

Suppose we want to test a coin for fairness. We flip the coin n times and record the number of times X0 that the coin turns up tails and the number of times X1 = n − X0 that the coin turns up heads. Now we set

$$

Z = ^1_{i=0}

$$

Then for a fair coin Z has approximately a chi-squared distribution with 2 − 1 = 1 degree of freedom. Verify this by computer simulation first for a fair coin (p = 1/2) and then for a biased coin (p = 1/3).

Answer:

First step: Simulating a coin toss with 1/2 chance of heads with 100 flips, doing this 1,000 times, and calculating x-squared based on the results of each simulation.

# Create a data frame to store the results of the simulations
results <- data.frame(matrix(ncol = 2, nrow = 1000))

# Create the first coin with 1/2 probability
coin_twosided <- c(0,1)

for (i in 1:1000) {
  # Flip coin 100 times and store a summary of the results in a table
  a <- table(sample(coin_twosided, size = 100, replace = TRUE))
  # Perform the chi square test so that we can get the X squared statistic (i.e. 'V')
  b <- chisq.test(a)
  # Save the statistic in a data frame
  results[i,1] <- b$statistic
}

Second step: Same as above but doing it with a “3-sided coin” (i.e. a 1/3 chance of getting any outcome)

# Create the first coin with 1/3 probability
coin_threesided <- c(0,1,2)

for (i in 1:1000) {
  # Flip coin 100 times and store a summary of the results in a table
  a <- table(sample(coin_threesided, size = 100, replace = TRUE))
  # Perform the chi square test so that we can get the X squared statistic (i.e. 'V')
  b <- chisq.test(a)
  # Save the statistic in a data frame
  results[i,2] <- b$statistic
}

Final step: Plot the results of the simulations for both coins and compare with chi square distributions. We can see that the distribution for each of the series of simulations is similar to the applicable chi squared distribution.

par(mfrow = c(2,2))
hist(results[,1], main="1/2 Probability Simulation")
curve(dchisq(x, df = 1), from = 0, to = 40)

hist(results[,2], main="1/3 Probability Simulation")
curve(dchisq(x, df = 2), from = 0, to = 40)