The Normal Distribution

Density (nv_dnorm), distribution function (nv_pnorm), quantile function (nv_qnorm), and random generation (nv_rnorm) for the Normal distribution with mean mean and standard deviation sd.

Usage

nv_dnorm(x, mean = 0, sd = 1, log = FALSE)

nv_pnorm(q, mean = 0, sd = 1, lower_tail = TRUE, log_p = FALSE)

nv_qnorm(p, mean = 0, sd = 1, lower_tail = TRUE, log_p = FALSE)

nv_rnorm(shape, mean = 0, sd = 1, dtype = NULL)

Arguments

x, q: (arrayish)
Quantiles at which to evaluate the density (x) or the distribution function (q).
mean: (arrayish)
Mean of the distribution. Either a scalar, or an array of exactly the shape of x/q/p (or the sample, for nv_rnorm).
sd: (arrayish)
Standard deviation of the distribution, shaped like mean. Must be positive, otherwise results are invalid.
log, log_p: (logical(1))
If TRUE, the densities/probabilities are given as logarithms. For nv_qnorm this describes the input p.
lower_tail: (logical(1))
If TRUE (default), probabilities are $P(X \le x)$; otherwise, $P(X > x)$.
p: (arrayish)
Probabilities at which to evaluate the quantile function. Values outside $[0, 1]$ give NaN.
shape: (integer())
Shape of the result.
dtype: (NULL | character(1) | DataType)
Floating point data type of the sample. The default (NULL) uses the common data type of mean and sd, and the default float type when both are R values.

Value

(arrayish)
nv_dnorm(), nv_pnorm() and nv_qnorm() return an arrayish with the shape of x/q/p and the float data type described under Data Types.

nv_rnorm() returns the sample, an arrayish of shape shape and the data type described under dtype.

Details

The Normal distribution has probability density function: $$f(x) = \frac{1}{\sigma\sqrt{2\pi}} \exp\left(-\frac{(x-\mu)^2}{2\sigma^2}\right)$$ where $\mu$ is the mean and $\sigma$ is the standard deviation.

nv_pnorm uses the asymptotic expansion from Abramowitz & Stegun (1964), equation 26.2.12, in the left tail when log_p = TRUE to maintain accuracy.

nv_qnorm uses the same minimax rational approximation as Moshier (1989) (this is ndtri in the Cephes library as used by JAX) for f64, and uses a new lower degree Remez minimax rational approximation on the same intervals for f32.

Data Types

nv_dnorm(), nv_pnorm() and nv_qnorm() compute at a float data type and return it. If x/q/p is a float, they compute at its data type; an R double settles on the default float. An integer x/q/p (an array or an R integer) is converted to the widest float among the arrays given for mean or sd, or to the default float if there is none. Integers beyond $2^{24}$ are not all exact at f32, so a large integer may change in that conversion.

An R value for mean or sd takes the computing data type. An array is promoted to it, so a wider one (e.g. f64 for an f32 x) is an error.

Random State

The draws come from the active random state: the global random state (see nv_set_seed()), or the random state that with_nv_rng_state() or local_nv_rng_state() make the default for a scope. They advance it in place.

Random generation

nv_rnorm samples via the Box-Muller transform. To sample with a covariance structure, use a Cholesky decomposition.

mean and sd are arrayish, so they may vary across the sample: they are applied to the draws after they have been reshaped to shape, and so may either be scalars or have exactly that shape. As in base R's rnorm(), an element whose mean is NaN, or whose sd is not finite or negative, is NaN. The random state is advanced regardless.

References

Abramowitz M, Stegun I (1964). Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, number 55 series Applied Mathematics Series. Dover Publications, New York. ISBN 0-486-61272-4.

Moshier S (1989). Methods and Programs for Mathematical Functions. Ellis Horwood. ISBN 0-7458-0289-3.

Examples

x <- nv_array(c(-1, 0, 1))
nv_dnorm(x)
#> AnvlArray
#>  0.2420
#>  0.3989
#>  0.2420
#> [ CPUf32{3} ] 
nv_dnorm(x, mean = 1, sd = 2)
#> AnvlArray
#>  0.1210
#>  0.1760
#>  0.1995
#> [ CPUf32{3} ] 
nv_dnorm(x, log = TRUE)
#> AnvlArray
#>  -1.4189
#>  -0.9189
#>  -1.4189
#> [ CPUf32{3} ] 

nv_pnorm(x)
#> AnvlArray
#>  0.1587
#>  0.5000
#>  0.8413
#> [ CPUf32{3} ] 
nv_pnorm(x, mean = 1, sd = 2)
#> AnvlArray
#>  0.1587
#>  0.3085
#>  0.5000
#> [ CPUf32{3} ] 
nv_pnorm(x, lower_tail = FALSE)
#> AnvlArray
#>  0.8413
#>  0.5000
#>  0.1587
#> [ CPUf32{3} ] 
nv_pnorm(x, log_p = TRUE)
#> AnvlArray
#>  -1.8410
#>  -0.6931
#>  -0.1728
#> [ CPUf32{3} ] 

p <- nv_array(c(0.025, 0.5, 0.975))
nv_qnorm(p)
#> AnvlArray
#>  -1.9600
#>   0.0000
#>   1.9600
#> [ CPUf32{3} ] 
nv_qnorm(p, mean = 1, sd = 2)
#> AnvlArray
#>  -2.9199
#>   1.0000
#>   4.9199
#> [ CPUf32{3} ] 
nv_qnorm(p, lower_tail = FALSE)
#> AnvlArray
#>   1.9600
#>  -0.0000
#>  -1.9600
#> [ CPUf32{3} ] 
nv_qnorm(nv_array(c(-700, -2, -0.1), dtype = "f64"), log_p = TRUE)
#> AnvlArray
#>  -37.2951
#>   -1.1015
#>    1.3096
#> [ CPUf64{3} ] 
with_nv_seed(42L, nv_rnorm(c(2, 3)))
#> AnvlArray
#>  -1.2172  1.9579  0.2044
#>  -0.8427 -0.1720  0.8960
#> [ CPUf32{2,3} ] 

# `sd` may also be an array of the same shape as the sample
sds <- nv_array(matrix(c(0.01, 0.1, 1, 10, 100, 1000), nrow = 2))
nv_rnorm(c(2, 3), sd = sds)
#> AnvlArray
#>    0.0087   1.0790  74.7185
#>    0.1132   5.1053 546.2689
#> [ CPUf32{2,3} ]