Tossing a coin in R

Imagine a coin. It has two sides. We assume that if we toss it it has exactly the same probability of falling on one or the other side. These sides are usually called Heads and Tails, but let’s simplify and use 1 and 0 instead.

We can create a coin object in R with only two values: 1 and 0

coin <- c(1,0)
coin
## [1] 1 0

How can we toss our coin? If we manually select one of the values, we would have to decide to select the first or the second, and that would do no good. Fortunately in R we have a sample function, that allows us to randomly select values from an object. Exactly what we need!

Remember that you can check the help page for this function with ?coin

sample(coin, size=1)
## [1] 0

Try that out a few times. Does it seem to be doing what we wanted? Remember that 1 represents Heads and 0 represents Tails

sample(coin, size=1)
## [1] 0
sample(coin, size=1)
## [1] 1
sample(coin, size=1)
## [1] 0

It’s much easier to change the size argument, which in this case simulates the number of times we are tossing our virtual coin. But, if we want to toss the coin multiple times we need to “reset” the coin each time, using the replace=TRUE argument. Otherwise our coin object will run out of values (by default sample will first randomly select one of the two values, the second time it will select the remaining value and it will not work a third time).

sample(coin, size=20, replace=TRUE)
##  [1] 1 1 1 1 1 1 0 0 0 1 1 1 0 0 1 0 0 0 1 0

So, this looks more like it. Although the above result is random (each time I run the code I get a different result), I’m pretty sure there are several 1 (Heads) and several 0 (Tails).

Exercises


Counting the number of Heads

Imagine you toss a coin 10 times.

How many Heads will you get? If you had to bet, what number would you choose? Maybe five?

tosses <- sample(coin, size=10, replace=TRUE)
sum(tosses)
## [1] 6

Oops, we got 6 Heads, you just lost your bet!

In the above example we decided to sum the results. Since we only have values of 1 and 0, summing essentially counts the number ones, which in our case are the Heads.

A more generic way would have been to say sum(tosses == 1), which would count the number of times the value 1 (in this case) appears in the object tosses. Do you remember the example in the first practical, where we did a sum(ages > 10)? If not, you might want to review that practical.

How come we got this value?

Remember we are simulating a random process. At each toss of a coin, we cannot know what side will fall. It could be Heads (1), it could be Tails (0), and each toss is independent from all other tosses.

If we use the code again, we will get a different value. Try it out, and you will see…

tosses <- sample(coin, size=10, replace=TRUE)
sum(tosses)
## [1] 2
tosses <- sample(coin, size=10, replace=TRUE)
sum(tosses)
## [1] 5
tosses <- sample(coin, size=10, replace=TRUE)
sum(tosses)
## [1] 6

There! We got a 5. But… also a 2 and a 6. How come?

Again, remember that we are simulating something random: the toss of a coin. Each of the 10 tosses can come up either 1 or 0 with equal probability. If you try to toss a real coin yourselves you would get something similar: perhaps 5 heads, but you could also get 4 or 6 or 3 or 7. It is in fact possible to get 0 or 10 as well, although you might think that would require a lot of “luck”. Also realize that the only impossible values are getting less than 0 or more than 10.

Exercises


Determining a probability

You might all agree that if you toss a coin, the probability of getting a Head is \(\frac{1}{2}\). But how do you come up with this idea?

We have seen in the previous exercises that it is not as easy to get a number of Heads that corresponds exactly to \(\frac{1}{2}\) of the number of times we tossed the coin. But this could simple be because we were using a small number of tosses.

In fact, one of the most common ways of defining probability consists of thinking of a process that you can repeat many, many, many times in exactly the same way.

Can you toss a coin many, many, many times? Of course.

Now, how many times should you get a Head if you tossed the coin a very very large number of times? As the number of tosses gets larger and larger, the fraction of Heads will get closer and closer to \(\frac{1}{2}\)

This concept is usually called the frequentist definition of probability. We can restate it as follows: the probability of an event is the fraction of times that event would occur if we were to repeat the same experiment an infinite number of times.

Now, we would never really be able to toss a coin an infinite number of times, but using R we can get a pretty good idea of what would happen, so let’s try it out!

Let’s toss our coin some more, how about 100 times:

tosses <- sample(coin, size=100, replace=TRUE)

And let’s calculate the fraction of Heads

sum(tosses) / 100
## [1] 0.48

Hmmm, pretty close, but not exactly 0.5. How about we try to toss the coin an even larger number of times?

tosses <- sample(coin, size=1000, replace=TRUE)
sum(tosses) / 1000
## [1] 0.472
tosses <- sample(coin, size=10000, replace=TRUE)
sum(tosses) / 10000
## [1] 0.5009
tosses <- sample(coin, size=100000, replace=TRUE)
sum(tosses) / 100000
## [1] 0.49926
tosses <- sample(coin, size=1000000, replace=TRUE)
sum(tosses) / 1000000
## [1] 0.500668

We are getting very close, but again not exactly 0.5, even when tossing our virtual coin one million times!

The values truly should be getting really close to 0.5. Perhaps a plot would be useful to see that this simulation process really does converge towards the true probability of getting a Head.

We need to again calculate different fractions of heads at different numbers of tosses, but this time saving the results for the plot.

I’ve tried to keep the code as simple as possible. Those of you who know a bit of R might be able to think of a better way to get the same result. For the rest of you: practice and you will get there…

tosses <- sample(coin, size=1000000, replace=TRUE)

frac1 <- sum(tosses[1])/1
frac2 <- sum(tosses[1:2])/2
frac3 <- sum(tosses[1:3])/3
frac4 <- sum(tosses[1:4])/4
frac5 <- sum(tosses[1:5])/5
frac6 <- sum(tosses[1:6])/6
frac7 <- sum(tosses[1:7])/7

frac_heads <- c(frac1, frac2, frac3, frac4, frac5, frac6, frac7)

A description of the above code: we first simulate tossing the coin 1 million times, and save the result. We then look at the first value, and calculate the fraction of heads. We then look at the first 2 values, and calculate the fraction of heads. We keep going until we’ve looked at the first 7 values (1:7). We then concatenate the 7 fractions and save them as frac_heads.

Finally, in the next bit of code we plot the seven fractions. The argument type allows to specify the type of plot (lines, points or both). The ylim argument specifies the limits of the y axis. The xlab and ylab arguments specify the labels of the x and y axes, respectively. The last line of code draws a horizontal line at 0.5, using the colour red.

plot(frac_heads, type="b", ylim=c(0,1), xlab="Number of tosses", ylab="Fraction of Heads")

abline(h=0.5, col='red')

This last plot shows what happens if you only toss the coin a few times. The results can easily differ quite a bit from 0.5.

To understand this behaviour, think of a few simple examples:

But we had already tossed our coin a million times. Let’s jump ahead and look at the fractions of Heads at much higher number of tosses

frac10 <- sum(tosses[1:10])/10
frac100 <- sum(tosses[1:100])/100
frac1k <- sum(tosses[1:1000])/1000
frac10k <- sum(tosses[1:10000])/10000
frac100k <- sum(tosses[1:1e5])/1e5
frac1000k <- sum(tosses[1:1e6])/1e6

frac_heads_big <- c(frac10, frac100, frac1k, frac10k, frac100k, frac1000k)

The code is almost the same as before (it’s a very good idea to reuse code!) but we use different names for all our objects. Note that we are still working with our same tosses object from above, where we had already saved 1 million tosses of our coin.

Note that it can be quite convenient to write certain numbers using scientific notation: 1e6 is exactly the same as \(10^6\) or 1 with 6 zeroes (1 million)

plot(frac_heads_big, type="b", ylim=c(0,1), xlab="Larger number of tosses", ylab="Fraction of Heads")

abline(h=0.5, col='red')

This last plot shows what happens when you toss the coin a much larger number of times. Now the fraction of heads is really getting close to the 0.5 value.

This result can also be used as an illustration of the Law of large numbers. This law states that with a large number of trials the average will be quite close to the expected value. The larger the number of trials, the closer it will be to the expected average. Realize that in our coin experiment, the expected average is also 0.5 (when half your numbers are 0 and half are 1, the average is 0.5).




Final Exercise


Imagine throwing dice. A die has six sides, instead of the two that a coin has. Let’s assume we are talking about a fair die, with exactly the same probability of falling on any side.

Let’s use R to create a virtual die and estimate the probability of “getting a 6”.

  • Create a virtual die in R, just like we did with a coin. Test it out by sampling one value multiple times.
  • Since we are only interested in the probability of getting a 6, can you think of different ways to create your virtual die to simplify this? Try it out…
  • Calculate and visualise the fraction of times your die lands on a 6 with more and more throws. What value does the probability converge towards?
  • How would you simulate the probability of getting “a 5 or a 6”?