Distributions and Random Numbers

For a random variable \( X \) with distribution function \( F \), the four functions answer four different questions, and the first argument of each one is named after the kind of number it takes:

Source: An Introduction to R: R as a set of statistical tables

The manual’s rule is to prefix the name of the family with d for the density, p for the cumulative distribution function, q for the quantile function and r for simulation, that is, random deviates. The first argument is x for the d functions, q for the p functions, p for the q functions and n for the r functions, with one exception family: rhyper(), rsignrank() and rwilcox() call the count of values nn, because their own parameter is already called n. The script runs the four normal functions side by side, shows the real argument names, and repeats the exercise for the binomial.

Source: R documentation: Normal

The binomial has a twist worth keeping in mind from the start: dbinom(3, size = 10, prob = 0.3) is the probability of exactly three successes (a number you can read as a probability), pbinom(3, ...) adds up the masses of 0, 1, 2 and 3, and qbinom(0.5, ...) returns an integer, the smallest count whose cumulative probability reaches one half. For a continuous distribution such as the normal, dnorm(0) is a density, 0.3989 here, and is not a probability: the probability of any single exact value is zero, and a density can be larger than one (the later script shows dnorm(0, sd = 0.1) is 3.99).

The stats package provides the families in the table below. Nineteen of them are the table in An Introduction to R; the multinomial, studentised range, birthday and Smirnov distributions come from the manual’s text and from the Distributions help page, which also lists the less common distributions of test statistics. The name of each family is what goes after the prefix, and the parameter names are the argument names after the first one (most functions also have the non-centrality parameter ncp where it makes sense; see the individual help pages).

Source: R documentation: Distributions

Not every family has all four prefixes. The multinomial distribution has dmultinom() and rmultinom() only; the studentised range (ptukey(), qtukey()) and the birthday problem (pbirthday(), qbirthday()) have only the p and q functions; and the two-sample Smirnov distribution has psmirnov(), qsmirnov() and rsmirnov() (new in R 4.2.0). The script checks every combination with exists() instead of trusting the table.

Source: R documentation: Distributions

Two things in the table deserve a sentence. The same family can have more than one parameterisation: the gamma takes rate or scale (the help page says scale defaults to 1 when it is omitted, and the two are related by scale = 1/rate), and the negative binomial takes prob or the mean mu (with prob = size/(size + mu)). And the geometric and negative binomial distributions count failures, not trials: dgeom(0, prob = 0.3) is 0.3, the probability of succeeding at the first try. The next script has one call per family, checks the two parameterisations against each other, and ends with the hypergeometric (drawing three books from a shelf of five poetry and eleven other books, the chance of at least one poetry book is 0.7054), the non-centrality parameter, and the distributions of test statistics.

Source: R documentation: NegBinomial

Distributions outside this list live in packages; the CRAN task view on probability distributions lists them (CRAN Task View: Probability Distributions), and the Distributions help page points to it.

The p and q functions all take two logical arguments, lower.tail and log.p, and the d functions take log. These are not cosmetic. lower.tail = FALSE asks for the upper tail \( P(X > q) \) directly, which is the quantity a one-sided p-value is; and the logarithmic versions keep numbers that would otherwise underflow or lose all their digits. The manual notes the cumulative hazard function as an example of what they are for: \( H(t) = -\log(1 - F(t)) \) is -pxxx(t, ..., lower.tail = FALSE, log.p = TRUE), and a log-likelihood is more accurate as the sum of dxxx(..., log = TRUE) than as the log of a product.

Source: An Introduction to R: R as a set of statistical tables

Every numeric argument except n is vectorised, and the result has the length of the longest one. That is more than a convenience: a vector of means or standard deviations gives a family of results in one call, which is how the tail probabilities of several models or the draws from several groups are computed without a loop. The rule is the ordinary recycling rule of R (shorter vectors are repeated; see Vectorisation and Iteration): pnorm(c(0, 1, 2, 3), mean = c(0, 10)) pairs the means 0, 10, 0, 10 with the four quantiles, and the second and fourth results are essentially zero because 1 and 3 are very far below a mean of 10. Only the first element of a logical argument is used, so lower.tail = c(FALSE, TRUE) means FALSE. For the r functions, n is the number of values to draw, but a vector of length greater than one is read as "this many values", length(n), so rnorm(c(10, 20, 30)) draws three numbers, which is a classic way to get a wrong sample size. Missing values propagate (pnorm(c(1, NA)) gives NA), and an NA parameter in an r function gives NaN with a warning, as shown in the errors script near the end.

R does not produce random numbers; it produces a deterministic sequence that behaves like random numbers. A generator is an algorithm with a state: each request for a random number reads the state, computes the next number and replaces the state. The sequence of numbers is therefore completely determined by the state at the start, and two sessions that start from the same state print the same "random" values, which is the property that makes simulations reproducible. Every distribution on this page is built from uniform numbers in the interval from 0 to 1 (the runif() stream): the normal generator by default transforms them by inversion, which is what qnorm() does (the help page for Random calls the default "Inversion"), the exponential, binomial and others use their own algorithms on top of the uniform stream, and sample() turns them into integers.

The state lives in an integer vector called .Random.seed in the global environment, and only there: the Random help page says the object is only looked for in the user’s workspace. With the default generator, the Mersenne-Twister, it has 626 integers: the first one codes the kinds of generators in force, the second is the current position in the state (it is 624 right after set.seed() and 1 after one draw) and the other 624 are the 32-bit words of the generator. The first element is a decimal code: the lowest two decimal digits give the uniform generator, the hundreds the normal generator and the ten thousands the sampler. In 10403 the 03 is the Mersenne-Twister (the codes are 0 for Wichmann-Hill, 1 for Marsaglia-Multicarry, 2 for Super-Duper, 3 for Mersenne-Twister, 4 for Knuth-TAOCP, 5 for a user-supplied generator, 6 for Knuth-TAOCP-2002 and 7 for L’Ecuyer-CMRG), the 4 is "Inversion" and the 1 is "Rejection". The script prints these facts and does not print the 624 words, which are not meant to be read. The help page says .Random.seed can be saved and restored but must not be edited by the user; restoring a saved copy, as the script does with assign(), is the supported use, and hand-editing the numbers is not.

Source: R documentation: Random

The object has no value before the first random number is needed. A fresh session does not hold a .Random.seed; the first call that needs a random number creates one from the current time and the process id, which is why two sessions that do not set a seed give different results. If a saved workspace is restored at startup, the seed is restored with it and the stream continues where the previous session stopped (see Sessions, Workspaces and Startup), another reason to set the seed explicitly at the top of a script that must be reproducible.

set.seed(seed, kind = NULL, normal.kind = NULL, sample.kind = NULL) is the recommended way to start a stream. It takes a single number, interpreted as an integer (a non-integer is truncated, so set.seed(123.9) is set.seed(123)), and expands it into the whole state of the generator in force, which makes it the simple way to get quite different states from small integers; the help page adds that there is no guarantee that different seeds give different streams, though exceptions would be extremely rare. The three other arguments choose the generators in the same call (see the next section), and the return value is NULL, invisibly. set.seed(NULL) re-initialises the generator as if no seed had ever been set, so the next draw is seeded from the clock again; it is the way to leave a reproducible stretch of a session and go back to unpredictable numbers.

Source: R documentation: Random

Two habits follow from this. Set the seed once, at the top of the script or at the start of each simulation, not inside a function that is called repeatedly (a set.seed(1) in every call makes every call return the same "random" numbers). And set it before the code whose output you want to reproduce, not after: the stream of a script depends on every earlier random draw, so a new rnorm(1) inserted above a plot changes everything below it. Seeding by hand from the clock, a habit found in old tutorials, is not needed: R already seeds from the clock when no seed was set, and set.seed(NULL) does it explicitly.

RNGkind(kind = NULL, normal.kind = NULL, sample.kind = NULL) queries or sets the three kinds in force, and returns the previous ones (invisibly when an argument is not NULL). NULL means "no change", and the word "default" returns to R’s default. The kinds are partially matched, so RNGkind("Mers") is the Mersenne-Twister. R 3.6.0 made the result three elements long, because the third one is the new sampler (before that, the result had two elements). The kinds of uniform generator are:

Source: R documentation: Random

The last column counts the integers of the seed; .Random.seed has one more, the leading code, and the script below prints the lengths. The default, the Mersenne-Twister, is what almost everybody should use; the help page lists the number of failures of each generator in the TestU01 Crush and BigCrush batteries (2 each for the Mersenne-Twister, 6 each for L’Ecuyer-CMRG, 40 clear failures of Crush for Marsaglia-Multicarry). The practical reasons to change the generator are reproducing an old result, and "L'Ecuyer-CMRG" for parallel work (see below). R 4.1.2 made selecting the Marsaglia-Multicarry generator warn in more places (PR#18168), because combined with some normal generators it shows deviations from normality.

The normal generator, normal.kind, is "Inversion" by default; the others are "Box-Muller" (stateful: it generates pairs and returns them one at a time, and the state is reset whenever the kind is selected), "Ahrens-Dieter", "Kinderman-Ramage", "Buggy Kinderman-Ramage" (the pre-1.7.0 version, with approximation errors, available only to reproduce old results and not accepted by set.seed()), and "user-supplied". The sampler, sample.kind, is "Rejection" by default since R 3.6.0 and "Rounding" before; see the section on sample(). The same seed with a different kind gives a different stream, as the script shows with one uniform draw per generator: the seed 42 gives 0.9148 under the Mersenne-Twister (the same number as set.seed(42); runif(1) gives) and 0.2508 under Wichmann-Hill. The script also lists the length of .Random.seed for each kind and the kind code in its first element.

Source: R documentation: Random

A user-supplied generator is a C function registered under the names that Random.user describes; it is the way to plug a hardware or library generator into R and is outside this page.

sample(x, size, replace = FALSE, prob = NULL) draws size elements from x, without replacement by default, so sample(x) is a random permutation of x. sample.int(n, size = n, replace = FALSE, prob = NULL, useHash = ...) is the bare version over the integers 1:n, where both n and size must be supplied as numbers and only the integers are drawn; it accepts an n larger than the biggest integer, up to the largest whole number a double can hold (the script draws from 1e10), and the result is then a double. The prob argument gives weights, which need not sum to one but must be non-negative and not all zero. With replacement, Walker’s alias method is used when there are more than 200 reasonably probable values; without replacement the weights are applied sequentially, the probability of the next element being proportional to the weights of those remaining, so the number of non-zero weights must be at least size. The useHash argument selects a hash-based algorithm that uses memory proportional to the sample and not to n; its default is n > 1e7 && !replace && is.null(prob) && size <= n/2. R 3.4.0 added useHash, and R 4.2.0 added sanity checks on size and n.

To sample rows of a data frame, sample their positions: books[sample(nrow(books), 3), ], or here the titles of three random books. To sample by group or with some rows more likely than others, use prob: the last block draws genres in proportion to their counts in books, 400 times with replacement, and the table of draws comes out close to the table of counts times 25 (that is 400/16: 50, 100, 50, 75, 50 and 75 expected, up to sampling noise).

Source: R documentation: sample

When x has length one, is numeric and is at least 1, sample(x) does not sample from x: it samples from 1:x. This convenience is documented, and the help page warns that this can give undesired behaviour when x has a varying length in calls such as sample(x). The bug appears when x is a computed subset that happens to have one element: sample(x[x > 9]) on x <- 1:10 is sample(10), a permutation of ten numbers, not a sample of the single element that was selected. The safe idiom is the one given in the help page’s examples: sample the positions with sample.int(), then index, x[sample.int(length(x), ...)]. Only a numeric vector of length one triggers the expansion; sample("only") returns "only", and non-integer numbers are truncated (sample(3.9) is a permutation of 1 to 3).

The errors of sample() have fixed messages, which the script collects: asking for more elements than there are without replacement is "cannot take a sample larger than the population when 'replace = FALSE'"; the wrong number of weights is "incorrect number of probabilities"; fewer positive weights than the sample size is "too few positive probabilities"; a negative weight is "negative probability"; and an empty or non-vector first argument is "invalid first argument".

Source: R documentation: sample

Before R 3.6.0, sample() turned a uniform number into an integer in 1:n by scaling it and rounding down, which makes sample() noticeably non-uniform on large populations; Ottoboni and Stark pointed this out, and the R NEWS records the change (PR#17494). From 3.6.0 the default sampler, "Rejection", draws random bits and rejects the values that would bias the result. The consequence is that the same seed gives different samples before and after: set.seed(1); sample(10) is 9 4 7 1 2 5 3 10 6 8 on R 3.6.0 or later and 3 4 5 7 2 8 9 6 10 1 on R 3.5.x or earlier, a difference visible in every tutorial or paper that quotes a sample drawn on R 3.5 or earlier. Only sample() and its relatives changed; the uniform and normal streams are the same (the script checks that runif(1) and rnorm(1) after the same seed are identical under both samplers).

To reproduce an old result, switch the sampler back: RNGkind(sample.kind = "Rounding"), or set.seed(seed, sample.kind = "Rounding"), or RNGversion("3.5.0"), which sets all three kinds as they were in that version. Each of them gives a warning, "non-uniform 'Rounding' sampler used", because the sampler is known to be biased; the script catches the warning to show its text. RNGversion() takes a version string and returns the same information as RNGkind() about the defaults of that version; for a version later than the running R it gives the current defaults, and RNGversion(getRversion()) works since R 3.6.0. For instance RNGversion("1.6.2") selects the Marsaglia-Multicarry generator with the Buggy Kinderman-Ramage normal generator, the defaults of R before 1.7.0, and RNGversion("1.7.0") gives the Mersenne-Twister with inversion and rounding. Restore the current generators afterwards with RNGkind(sample.kind = "default") or RNGversion("4.6.1"): leaving "Rounding" on silently biases every later sample() in the session.

Source: R NEWS: Changes in R 3.6.0

R-devel, the development version that will become R 4.7.0, has a change to the random binomial generator that affects reproducibility; it is described here as upcoming and none of it is run on this page. The R-devel NEWS says that the BTPE algorithm of rbinom() is fixed, following a report from Kate Hudson (PR#19049), and that to keep old results reproducible RNGkind() gets a new argument binom.kind, with RNGversion("4.6") keeping the previous behaviour. If it ships as announced, a script that draws binomial numbers under a fixed seed on R 4.7.0 may print different numbers from R 4.6.1, and the way to pin the old ones is the same as for sampling in 3.6.0: RNGversion("4.6") at the top of the script. The details may change before release; the NEWS of the release is the reference.

Source: R NEWS: Changes in R-devel

Simulation means drawing random data, computing a statistic and repeating many times. replicate(n, expr, simplify = "array") is the base R tool: the help page describes it as a wrapper for the common use of sapply() to evaluate an expression repeatedly, and its body is a single call, which the first lines of the script print. The key point is that expr is an expression, re-evaluated on every repetition, not a value computed once; replicate(5, runif(1)) gives five different numbers, while rep(runif(1), 5) repeats one. A single value per repetition gives a vector, a vector of length k gives a matrix with k rows and n columns, simplify = FALSE gives a list and simplify = "array" (the default) can build higher-rank arrays from matrices. The strict alternative is vapply() over seq_len(n), which checks the type and length of every result; with the same seed it draws the same numbers as replicate() (the script checks that the two vectors are identical), so the choice is a matter of clarity and strictness, not of speed or of the stream.

Source: R documentation: lapply, sapply, vapply and replicate

Monte Carlo estimation replaces a probability that is hard to compute by the share of simulated experiments in which the event happens. The script estimates the probability that a fair die shows a six at least twice in ten rolls from 20000 simulated experiments: 0.5134, against the exact value 0.5155 from pbinom(). The error of the estimate is about the square root of \( p(1-p)/N \), here about 0.0035, so agreement to about two decimals is what to expect; to halve the error, quadruple the number of repetitions. The exact value is available here because the experiment is a binomial; for an event that no formula covers, the simulation is the method. The pattern to copy is the one-line experiment: sample(6, 10, replace = TRUE) rolls ten dice at once, the comparison and sum() count the sixes, and mean() of the logical results is the estimated probability.

The bootstrap estimates the sampling distribution of a statistic from the sample itself: draw a new sample of the same size with replacement from the data, compute the statistic, repeat thousands of times, and look at the spread of the results. The shelf example bootstraps the mean of books$pages: 14 of the 16 books have a known page count, their mean is 268.2 pages, and 5000 resamples give bootstrap means centred at 268.1 with a standard deviation of 27.8 (the textbook standard error sd(pages) / sqrt(n) is 28.9) and a percentile interval of 213.6 to 322.4 pages, the 2.5th and 97.5th percentiles of the bootstrap means. Two details of the script matter. The vapply form with sample.int(length(pages), replace = TRUE) draws exactly the same numbers as the replicate() form (the script prints TRUE), and it is the safe one if pages could ever shrink to one element, because of the length-one trap above. And replicate() can return several statistics per repetition, a matrix with one row per statistic, so the mean and the median of each resample come out of one pass. With 14 observations the interval is wide and the percentile interval is rough, a reminder that the bootstrap does not create information the sample does not have.

Source: R documentation: sample

Counts of events are modelled by the Poisson distribution when the events are independent and have a constant rate, and by the binomial when there is a fixed number of trials. The shelf has 59 loans with a valid date, spread over the 72 weeks counted from the first loan to the last, and the number of loans per week has mean 0.82 and variance 0.80: the mean and the variance of a Poisson variable are equal, and these are almost equal, a first sign that the Poisson model is reasonable. The expected number of weeks with 0, 1, 2, 3 and 4 loans under a Poisson model with that mean, length(per_week) * dpois(0:4, lambda), is 31.7, 26.0, 10.7, 2.9 and 0.6 weeks, against the observed 31, 27, 11, 2 and 1: a good fit. A simulated year of the same shelf is rpois(72, lambda = mean(per_week)). The other generators of the script are the ones used the most: rbinom() for a number of successes in a fixed number of trials; rexp() for waiting times (a mean gap of 7 days between loans is a rate of 1/7); runif() for a number in an interval; rnorm() with a vector of means, which draws one value around each mean (three means give three draws, not three hundred; the second argument may also be a vector of standard deviations); and rmultinom(), which returns a matrix with one column per draw and one row per category, here three draws of ten loans spread over three genres, with dmultinom() giving the probability of one particular split (0.085 for 5, 3, 2).

Source: R documentation: Multinom

The integer-valued generators rbinom(), rgeom(), rhyper(), rpois(), rnbinom(), rsignrank() and rwilcox() return integer vectors, as they have since R 3.0.0, and since R 4.0.0 they return doubles, with no NA, when the numbers would be outside the integer range.

A permutation test asks how unusual the observed difference between two groups is if the group labels were irrelevant: shuffle the labels with sample(), recompute the difference, repeat, and see where the observed value falls. The shelf example compares the loan durations of long books (300 pages or more: 17 of the 42 returned loans of books with a known page count) with the others. The observed difference of means is 0.22 days, 94 per cent of the 10000 shuffles give a difference at least as large in size (a two-sided permutation p-value of 0.9435), and the classical Welch t test gives 0.9342, the same conclusion by a formula: there is no evidence that the length of a book changes how long it is kept. The permutation test needs no distribution assumption, only the exchangeability of labels under the null, and it uses nothing but sample() on the labels; the tests themselves, and their assumptions, are in Descriptive Statistics and Tests.

Source: R documentation: t.test