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.
