shotGroups (version 0.7.1)

mvnEll: 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)
qmvnEll(p,   sigma = diag(2), mu, e, x0, lower.tail = TRUE, loUp=NULL)
rmvnEll(n,   sigma = diag(2), mu, e, x0, method = c('eigen', 'chol', 'cdf'), loUp=NULL)

Arguments

r

vector of radii for the offset ellipse.

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 matrix 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]$.

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.

Details

pmvnEll is implemented by first transforming the integration region to the unit disc/sphere, then decorrelating the normal distribution through rotation, and finally applying farebrother to calculate the quadratic form.

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.

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.

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"))