Vectorised statistics on distributions commonly require two inputs:
the distribution(s), and the point(s) at which the statistic is
evaluated. The distributional package offers several broadcasting
mechanisms for statistical operations on distributions, which are safer
alternatives to the recycling rules of base R’s
p/d/q/r distribution
functions.
dist <- dist_normal(mu = c(0, 3), sigma = c(1, 2))
dist
#> <distribution[2]>
#> [1] N(0, 1) N(3, 4)The form of the argument given to the statistic
(e.g. density(), cdf(),
quantile(), or hilo(), or
generate()) changes how the statistic is applied to the
distributions.
| Argument | Combined by | Result |
|---|---|---|
| Scalar | Broadcasting to every distribution | Unwrapped, one result per distribution |
| Vector | Mapping every point onto every distribution |
list, one element per distribution |
| List | Recycling each element against the distributions |
data.frame, one column per element |
density(dist, at = 0) # broadcast
#> [1] 0.3989423 0.0647588
density(dist, at = c(0, 1)) # mapped
#> [[1]]
#> [1] 0.3989423 0.2419707
#>
#> [[2]]
#> [1] 0.0647588 0.1209854
density(dist, at = list(d = c(0, 1))) # recycled
#> d
#> 1 0.3989423
#> 2 0.1209854In most cases, a scalar or vector is the natural choice - the same
statistic is computed for every distribution. List arguments are useful
when you want to vary the statistic’s arguments across distributions,
allowing for similar behaviour to the
p/d/q/r distribution
functions of base R.
A single point is evaluated for every distribution, giving one value each. The result is an atomic vector with the same length as the distributions:
density(dist, 0)
#> [1] 0.3989423 0.0647588A bare vector is a set of points, and every point is evaluated against every distribution. Two distributions and two points give all four values, as one element per distribution:
The first element holds the density of N(0, 1) at
0 and 1, the second the density of
N(3, 4) at the same two points.
Because each distribution is handled separately, the number of points is unrelated to the number of distributions. This makes a dense grid against a single distribution natural, as when drawing a density curve or a quantile dotplot:
quantile(dist_normal(0, 1), ppoints(5))
#> [[1]]
#> [1] -1.1797611 -0.4972006 0.0000000 0.4972006 1.1797611Broadcasting a scalar is the degenerate case of this: one point per distribution, with the list of length-1 vectors simplified to a bare vector.
A list argument for statistics behaves differently, each list element
describes statistics computed from each distribution. This is closer to
the behaviour of base R’s
p/d/q/r functions,
where different statistical arguments apply to different
distributions.
Consider dnorm() for the density of two normal
distributions at two points:
This computes the density of N(0, 1) at 0
and of N(3, 4) at 1. Almost always it is a
mistake to compute densities at different points on each distribution
(the mapping behaviour above is usually needed), but for full control
over the arguments, a list can be used to specify the arguments for each
distribution:
These are the two values on the diagonal of the mapped result above:
the first distribution is evaluated only at 0, the second
only at 1.
The result is a data frame with one row per distribution and one column per list element. The list elements are recycled against the distributions, so a single value is used for every distribution:
cdf(dist, list(fixed = 0, varying = c(0, 1)))
#> fixed varying
#> 1 0.5000000 0.5000000
#> 2 0.0668072 0.1586553The recycling behaviour follows the vctrs recycling rules, where only size 1 arguments are recycled to match the number of distributions.
quantile(dist, list(p = c(0.1, 0.5, 0.9)))
#> Error in `FUN()`:
#> ! Cannot recycle input of size 3 to match the distributions (size 2).This is safer than the partial recycling behaviour that the
p/d/q/r functions in
base R use, which is a common source of silent errors in analysis.
Three points against two distributions silently wraps around to give
three values: N(0, 1) at -1,
N(3, 4) at 0, and then N(0, 1)
again at 2.
The same rules apply to multivariate distributions, but with matrices rather than vectors. Each column of the matrix relates to a variate of the distributions. This makes a 1-row matrix the “scalar” and a multi-row matrix the “vector” of the rules above.
mvn <- dist_multivariate_normal(
mu = list(c(0, 0), c(3, 3)),
sigma = list(diag(2), 2 * diag(2))
)
density(mvn, cbind(0, 0)) # 1-row matrix: broadcast
#> [1] 0.1591549431 0.0008840259
density(mvn, cbind(c(0, 1), c(0, 1))) # 2-row matrix: mapped
#> [[1]]
#> [1] 0.15915494 0.05854983
#>
#> [[2]]
#> [1] 0.0008840259 0.0107696397Lists of matrices are recycled against the distributions.