CODEWITHSUNDEEP
rloops
Antoni Parellada
Antoni Parellada

Reputation: 4791

Simulating coin toss

In the New York Times yesterday there was a reference to a paper essentially saying that the probability of 'heads' after a 'head' appears is not 0.5 (assuming a fair coin), challenging the "hot hand" myth. I want to prove it to myself.

Thus, I am working on coding a simulation of 7 coin tosses, and counting the number of heads after the first head, provided, naturally, that there is a first head at all.

I came up with the following lines of R code, but I'm still getting NA values, and would appreciate some help:

n <- 7              # number of tosses
p <- 0.5            # probability of heads
sims <- 100         # number of simulations

Freq_post_H <- 0    # frequency of 'head'-s after first 'head' 
    for(i in 1:sims){
        z <- rbinom(n, 1, p)
        if(sum(z==1)!=0){
        y <- which(z==1)[1]
        Freq_post_H[i] <- sum(z[(y+1):n])/length((y+1):n) 
        }else{
            next()
        }
    Freq_post_H
    }
Freq_post_H

What am I missing?

CONCLUSION: After the initial hiccups of mismatched variable names, both responses solve the question. One of the answers corrects problems in the initial code related to what happens with the last toss (i + 1) by introducing min(y + 1, n), and corrects the basic misunderstanding of next within a loop generating NA for skipped iterations. So thank you (+1).

Critically, and the reason for this appended "conclusion" the second response addresses a more fundamental or conceptual problem: we want to calculate the fraction of H's that are preceded by a H, as opposed to p(H) in whatever number of tosses remain after a head has appeared, which will be 0.5 for a fair coin.

Upvotes: 2

Views: 9386

Answers (4)

Jon Dollard
Jon Dollard

Reputation: 1

Below is some sample code in R to simulate a fair coin toss in R using the sample function. You can modify it as you like to simulate any number of flips. Since the outcome of flipping a coin is independent for each flip, the probability of a head or tail is always 0.5 for any given flip. Over many coin flips the probability of at least half of the flips being heads (or tails) will converge to 0.5. The probability that you get exactly half heads and half tails approaches 0.

n <- 7
count_heads <- 0
coin_flip <- sample(c(0,1), n, replace = TRUE)
for(flip_i in 1:n)
{
  if(coin_flip[flip_i] == 1)
  {
    count_heads = count_heads + 1
  }
}
count_heads/n

Upvotes: 0

SATHPra
SATHPra

Reputation: 1

This is to simulate 100 fair coin tosses 30,000 times

counter <- 1 
coin <- sum(rbinom(100,1,0.5)) 
while(counter<30000){
  coin <- c(coin, sum(rbinom(100,1,0.5)))
  counter <- counter+1
}

Try these after running above variable

 hist(coin)
 str(coin)
 mean(coin)
 sd(coin)

Upvotes: 0

Roland
Roland

Reputation: 132706

This is a simulation of what they did in the newspaper:

nsims <- 10000
k <- 4
set.seed(42)
sims <- replicate(nsims, {
  x <- sample(0:1, k, TRUE)
  #print(x)
  sum( # sum logical values, i.e. 0/1
   diff(x) == 0L & # is difference between consecutive values 0? 
     x[-1] == 1L ) / # and are these values heads? 
       sum(head(x, -1) == 1L) #divide by number of heads (without last toss)
})

mean(sims, na.rm = TRUE)  #NaN cases are samples without heads, i.e. 0/0
#[1] 0.4054715

k <- 7

sims <- replicate(nsims, {
  x <- sample(0:1, k, TRUE)
  #print(x)
  sum(diff(x) == 0L & x[-1] == 1L) / sum(head(x, -1) == 1L) 
})

mean(sims, na.rm = TRUE) 
#[1] 0.4289402

Upvotes: 3

Pierre L
Pierre L

Reputation: 28441

n <- 7              # number of tosses
p <- 0.5            # probability of heads
sims <- 100         # number of simulations

Prob_post_H <- 0    # frequency of 'head'-s after first 'head' 
    for(i in 1:sims){
        z <- rbinom(n, 1, p)
        if(sum(z==1) != 0){
        y <- which(z==1)[1]
        Prob_post_H[i] <- mean(z[min(y+1, n):n], na.rm=TRUE)
        }else{
            next()
        }
    }
mean(Prob_post_H,na.rm=TRUE)
#[1] 0.495068

It looks like it's right around 50%. We can scale up to see more simulations.

sims <- 10000

mean(Prob_post_H,na.rm=TRUE)
#[1] 0.5057866

Still around 50%.

Upvotes: 1

Related Questions