# Simulating the outcome of X number of random trials

Let's say I want to flip a coin a billion times and count how many heads I get. But I don't have time to actually flip a coin a billion times.

Assuming I have a random number generator than can generate a number between 0 and 1, is there a technique I can use to simulate the outcome of 1 billion random coin tosses, such that it's statistically equivalent to actually doing a billion random coin tosses?

I have been thinking about this a lot and Googling, but haven't come up with a solution. It feels to me like we could plot a curve like this (please forgive my horrible sketch): Where x is the chance of getting fewer than f(x) heads.

So when x is 0 (0% chance), f(x) is also 0. There's a 0% chance of getting fewer than 0 heads.

When x is 0.5 (50% chance), f(x) is half a billion. There's a 50% chance of getting fewer than half a billion heads.

When x is 1 (100% chance), f(x) is a billion (or a billion and one, I suppose). There's a 100% chance of getting fewer than a billion heads.

Moving right from 0 or left from 1, f(x) very quickly goes to very close to half a billion. This reflects how it would be extraordinarily unlikely to get, say, 400 million or 600 million heads.

My thinking is that if I had a correct equation for f(x), then I could generate a random number between 0 and 1, calculate f(x), and that would be the number of heads from my one billion "trials".

Would this concept actually work and is there a known approach (even if completely different from what I'm thinking) for solving this problem?

What you apparently have in mind is known as the inverse transform sampling method. This is based on the cumulative distribution function (cdf) of the distribution you want to generate. This is what you call f, and in a graph as a standard the potential outcomes would be on the x- and the cdf values on the y-axis. You did it the other way round, which isn't much of a problem; I just say this because you may want to compare this with Binomial cdf graphs you may find somewhere. The distribution you need is a Binomial with parameters "one Billion" (number of trials) and 0.5 (success probability for each trial). I haven't thought about whether you will run into numerical problems evaluating this for $$n=10^9$$; in that case there should be approximations.

• R falls back to inverse sampling when the number of trials is large, so I guess it should work at least to some degree github.com/wch/r-source/blob/…
– Tim
Jun 4 at 18:26
• At stats.stackexchange.com/a/608790/919 I show how to generate samples of size $10^{30}$ or larger and at stats.stackexchange.com/a/606853/919, a sample of size $10^{100}$ is shown. The ideas applied to these more complicated problems translate directly to this Binomial sampling problem.
– whuber
Jun 5 at 13:15
• This does look like what I need. Appreciated! Jun 5 at 21:11

The distribution is Binomial, which is discrete. Using the generic function s given at What is the statistical distribution of the number of times a number i would appear after M trials of a poisson distribution? for drawing random values from any discrete distribution, we can sample such a large number of coin flips extremely efficiently.

As an example, the implementation rBinomial (given below) will repeat this experiment of flipping a coin one billion times as often as you like. For instance, the call rBinomial(1e20, 1e9) repeats this experiment $$10^{20}$$ times but takes only two seconds to execute. Because it returns almost a quarter million distinct values, which can take a long time to plot, I have plotted a random sample of a thousand of these values showing how many times each value was the result of the experiment: Of course this is identical to the rounded values of a Normal distribution: but the point is that the exact probabilities (in double precision floating point, of course) were used. Moreover, if you tried to use the Normal approximation in this instance it would fail, because the software simply cannot generate $$10^{20}$$ values in reasonable time.

The short answer to the question of flipping a coin one billion times, then, is obtained by computing

rBinomial(1, 1e9)


This takes a half second (not bad for a billion flips!) and returns a random value such as

500015061

The arguments to rBinomial are

• n is the number of times to replicate the coin-flipping experiment.
• N is the number of times to flip the coin in each experiment.
• p is the chance of heads on each flip.
• Z.max helps limit the storage used by the algorithm. It is a Z-score. By default it's effectively no limit at all (because the binomial probability associated with larger Z-scores is vanishingly small).
• k.max determines when to resort to a Normal approximation. When the range of likely values exceeds k.max, the values will be found by rounding Normal variates. This will work only when n is sufficiently small (less than $$10^8$$ or so is practicable).

Its output is a two-column array of values and their counts. All counts are positive and (therefore) the values are not necessarily consecutive. For example, one possible output of rBinomial(5, 10) is

     Value Count
[1,]     2     1
[2,]     5     2
[3,]     6     2


This documents five repetitions of ten coin flips. In one case there were two heads, in two cases there were five heads, and in two cases there were six heads.

rBinomial <- function(n, N, p = 1/2, Z.max = 8, k.max = 1e6) {
Z <- min(Z.max, -qnorm(1e-3 / N)) # (No need to look beyond this z-score)
mu <- p * N
sigma <- exp((log(p) + log(1 - p) + log(N)) / 2)
x.0 <- max(0, floor(mu - Z*sigma))
x.1 <- min(N, ceiling(mu + Z*sigma))
if (x.1 - x.0 >= k.max) {
warning("Normal approximation employed.")
round(rnorm(n, mu, sigma)) # Use Normal approximation
} else {
x <- seq(x.0, x.1, by = 1)
probs <- dbinom(x, N, p)
X <- cbind(Value = x, Count = s(n, probs))
i <- X[, 2] > 0
X[i, , drop = FALSE]
}
}
#
# Flip a coin one billion (1e9) times.
#
system.time(x <- rBinomial(1, 1e9))
#
# Repeat this process 1e20 (100,000,000,000,000,000,000) times.
#
system.time(x <- rBinomial(1e20, 1e9)) # Around 2 seconds on a 2021 workstation
#
# Sketch the results.  The output is an array with 'Value' and 'Count' columns
# tallying all 1e20 experimental results.
#
j <- sample.int(nrow(x), 1e3)
plot(x[j, ], type = "h")


The function s referenced in this code can be found at https://stats.stackexchange.com/a/606853/919.

Define a random variable $$X$$, such that $$X=0$$ if the coin is landing with a head, or $$X=1$$ if the coin is landing with a tail. Suppose you want to simulate the cases with probability $$p$$ for $$X=0$$.

Use a random number generator to generate i.i.d. uniformly distributed numbers between $$0$$ and $$1$$, denoted by $$S$$.

Then, quantize $$S$$ to $$X$$, that is, corresponding to each $$S$$, output $$X=0$$ if $$S, and $$X=1$$ if $$S>p$$.

In rare case if $$S=p$$ exactly, simply convert $$S$$ to $$X=0$$ and $$X=1$$ alternately for fairness.

Then, each value of $$X$$ represents the result of flipping a coin.

• The point of the question is to generate huge samples. Replacing one "billion" with $10^{30},$ for instance, shows this approach won't work.
– whuber
Jun 5 at 13:16