F_from_f {sphunif}

Distribution and quantile functions from angular function

Description

Numerical computation of the distribution function F and the quantile function F^{-1} for an angular function f in a tangent-normal decomposition. F^{-1}(x) results from the inversion of

F(x) = \int_{-1}^x \omega_{p - 1}c_f f(z) (1 - z^2)^{(p - 3) / 2} \,\mathrm{d}z

for x\in [-1, 1], where c_f is a normalizing constant and \omega_{p - 1} is the surface area of S^{p - 2}.

Usage

F_from_f(f, p, Gauss = TRUE, N = 320, K = 1000, tol = 1e-06, ...)

F_inv_from_f(f, p, Gauss = TRUE, N = 320, K = 1000, tol = 1e-06, ...)

Arguments

f

angular function defined on [-1, 1]. Must be vectorized.

p

integer giving the dimension of the ambient space R^p that contains S^{p-1}.

Gauss

use a Gauss–Legendre quadrature rule to integrate f with N nodes? Otherwise, rely on integrate Defaults to TRUE.

N

number of points used in the Gauss–Legendre quadrature. Defaults to 320.

K

number of equispaced points on [-1, 1] used for evaluating F^{-1} and then interpolating. Defaults to 1e3.

tol

tolerance passed to uniroot for the inversion of F. Also, passed to integrate's rel.tol and abs.tol if Gauss = FALSE. Defaults to 1e-6.

...

further parameters passed to f.

Details

The normalizing constant c_f is such that F(1) = 1. It does not need to be part of f as it is computed internally.

Interpolation is performed by a monotone cubic spline. Gauss = TRUE yields more accurate results, at expenses of a heavier computation.

If f yields negative values, these are silently truncated to zero.

Value

A splinefun object ready to evaluate F or F^{-1}, as specified.

Examples

f <- function(x) rep(1, length(x))
plot(F_from_f(f = f, p = 4, Gauss = TRUE), ylab = "F(x)", xlim = c(-1, 1))
plot(F_from_f(f = f, p = 4, Gauss = FALSE), col = 2, add = TRUE,
     xlim = c(-1, 1))
curve(p_proj_unif(x = x, p = 4), col = 3, add = TRUE, n = 300)
plot(F_inv_from_f(f = f, p = 4, Gauss = TRUE), ylab = "F^{-1}(x)")
plot(F_inv_from_f(f = f, p = 4, Gauss = FALSE), col = 2, add = TRUE)
curve(q_proj_unif(u = x, p = 4), col = 3, add = TRUE, n = 300)