mvnEll: Multivariate normal offset ellipse probabilities

In dwoll/shotGroups: Analyze Shot Group Data

Multivariate normal offset ellipse probabilities

Description

Probability of an offset ellipsoid for a correlated multivariate normal distribution. Offset circle probabilities are a special case.

Usage

pmvnEll(r=1, sigma = diag(2), mu, e, x0, lower.tail = TRUE,
        method_cdf = c('integrate', 'saddlepoint'))
        
qmvnEll(p,   sigma = diag(2), mu, e, x0, lower.tail = TRUE,
        loUp=NULL, method_cdf = c('integrate', 'saddlepoint'))

rmvnEll(n,   sigma = diag(2), mu, e, x0,
        method = c('eigen', 'chol', 'cdf'), loUp=NULL,
        method_cdf = c('integrate', 'saddlepoint'))

Arguments

r	vector of radii for the offset ellipse defined by e.
p	vector of probabilities.
n	number of observations. If length(n) > 1, the length is taken to be the number required.
sigma	true positive definite covariance matrix of multivariate normal distribution.
mu	true center of multivariate normal distribution.
e	positive definite matrix characterizing the offset ellipse defined by (x-x0)' e (x-x0) < r^2. If the ellipse defined by e has semi-axis lengths equal to the square root of the eigenvalues of a matrix S, and is oriented along the eigenvectors of S, then e = S^-1. By default a circle.
x0	center of the offset ellipse.
method	string indicating which method to use for generating random deviates. See details.
loUp	search interval for numerical root finding. Either a vector with the lower and upper interval boundary, a list of such vectors, or an (n x 2)-matrix. See details.
lower.tail	logical. If TRUE (default), probabilities are P[X \le x] otherwise, P[X > x].
method_cdf	string indicating which method to use for calculating the sum of non-central chi^2 variables. See details.

Details

pmvnEll is implemented by first transforming the integration region to the unit disc/sphere, then decorrelating the normal distribution through rotation. Finally, the quadratic form (sum of non-central chi^2 variables) is calculated depending on method_cdf. For method_cdf='integrate', numerical integration using farebrother is chosen, for method_cdf='saddlepoint', Kuonen's (1999) saddlepoint approximation. Note that the equation for the cumulant generating function K has a missing zeta in the numerator of the second term in Kuonen (1999), see Imhof (1961) instead. Lower tail probabilities are calculated as '1 - upper tail probability', so loss of accuracy is likely when the upper tail probability is very small (<1e-9).

qmvnEll is implemented through numerical root finding of pmvnEll. If no search interval for uniroot is provided, the quantiles of an approximating non-central chi^2 distribution are used to determine the search intervals.

rmvnEll with method='eigen' or with method='chol' simulates 2D normal deviates based on sigma and mu, and then determines the radius around x0. rmvnEll with method='cdf' is much slower as it performs numerical root finding of pmvnEll given simulated quantiles from a uniform random variable in (0,1). If no search interval for uniroot is provided, the quantiles of an approximating non-central chi^2 distribution are used to determine the search intervals.

See Hoyt for the distribution of radial error around the true center of correlated bivariate normal variables with unequal variances. See Rice for the distribution of radial error around an offset center for uncorrelated bivariate normal variables with equal variances. See Rayleigh for the distribution of radial error around the true center of uncorrelated bivariate normal variables with equal variances.

Value

pmvnEll integrates the multivariate normal distribution over an arbitrary ellipsoid and thus gives the cumulative distribution function. qmvnEll gives the quantile function, rmvnEll generates random deviates.

The functions are vectorized in r and p but not in the remaining parameters.

References

DiDonato, A. R., & Jarnagin, M. P. (1961a). Integration of the general bivariate Gaussian distribution over an offset circle. Mathematics of Computation, 15 (76), 375-382.

DiDonato, A. R., & Jarnagin, M. P. (1961b). Integration of the general bivariate Gaussian distribution over an offset ellipse (NWL TR 1710). Dahlgren, VA: U.S. Naval Weapons Laboratory.

Duchesne, P., & Lafaye de Micheaux, P. (2010). Computing the distribution of quadratic forms: Further comparisons between the Liu-Tang-Zhang approximation and exact methods. Computational Statistics and Data Analysis, 54, 858-862.

Imhof, J. P. (1961). Computing the distribution of quadratic forms in normal variables. Biometrika, 48, 419-426.

Kuonen D. (1999). Saddlepoint Approximations for Distributions of Quadratic Forms in Normal Variables. Biometrika, 86, 929-935.

See Also

Hoyt, farebrother, uniroot

Examples

# define a bivariate normal distribution
mu    <- c(2, -1)                        # true mean
sigma <- cbind(c(10, 6), c(6, 10))       # covariance matrix

# define circular integration region
ctr <- c(1, 0)                           # center
e1  <- diag(2)                           # circle
r   <- 2                                 # radius
pmvnEll(r,   sigma=sigma, mu=mu, e=e1, x0=ctr) # probability
qmvnEll(0.5, sigma=sigma, mu=mu, e=e1, x0=ctr) # quantile
rmvnEll(5,   sigma=sigma, mu=mu, e=e1, x0=ctr) # random numbers

# define elliptical integration region
S  <- cbind(c(3.5, -0.3), c(-0.3, 1.7))
e2 <- solve(S)
pmvnEll(r,   sigma=sigma, mu=mu, e=e2, x0=ctr) # probability
qmvnEll(0.5, sigma=sigma, mu=mu, e=e2, x0=ctr) # quantile
rmvnEll(5,   sigma=sigma, mu=mu, e=e2, x0=ctr) # random numbers

# plot all regions
evSig <- eigen(sigma)$values
evS   <- eigen(S)$values
xLims <- range(c( mu[1]+c(-1, 1)*sqrt(evSig[1]),
                 ctr[1]+c(-r, r)*sqrt(evS[1])))
yLims <- range(c( mu[2]+c(-1.25, 1.25)*sqrt(evSig[1]),
                 ctr[2]+c(-r, r)*sqrt(evS[1])))

plot(xLims, yLims, type="n", asp=1)
points(mu[1],  mu[2],  pch=16, cex=2, col="black")
points(ctr[1], ctr[2], pch=15, cex=2, col="blue")
drawEllipse(mu, sigma, r=0.75, fg="black")
drawEllipse(mu, sigma, r=1,    fg="black")
drawEllipse(mu, sigma, r=1.25, fg="black")
drawEllipse(mu, sigma, r=1.5,  fg="black")
drawEllipse(ctr, e1, r=r, fg="blue")
drawEllipse(ctr, S,  r=r, fg="red")
legend(x="bottomright", legend=c("normal iso-densities",
       "integration circle", "integration ellipse"),
       lty=1, col=c("black", "blue", "red"))

dwoll/shotGroups documentation built on Feb. 16, 2024, 2:21 p.m.