13
$\begingroup$

In a real-valued multivariate case, is there a way to uniformly sample the points from the surface where the Mahalanobis distance from the mean of the is a constant?

EDIT: This just boils down to sampling points uniformly from the surface of a hyper-ellipsoid that satisfies the equation,

$$(x-\mu)^T \Sigma^{-1}(x-\mu) = d^2.$$

To be more precise, by "uniformly", I mean sample such that each area element $dA$ of the hyper-surface contains the same probability mass.

$\endgroup$
16
  • 1
    $\begingroup$ Correct me if I'm wrong: are you asking "given a random variable $ X $, how can I uniformly sample from the points that are a given Mahalanobis distance $ c $ away from $ \mathbb{E}[X] $?" $\endgroup$
    – Kevin Li
    Sep 16, 2018 at 22:17
  • 1
    $\begingroup$ I think we will need a suitable definition of "uniformly." The reason is this: in two dimensions, this set of points lies along some ellipse. Is one supposed to sample from that ellipse in such a way that equal lengths have equal chances, or that equal angles have equal chances, or so that equal lengths when the variables are standardized have equal chances, or in some other way? If you could explain what this sampling aims to achieve, that might give us enough information to know what you are trying to ask. $\endgroup$
    – whuber
    Sep 17, 2018 at 14:07
  • 1
    $\begingroup$ These surfaces are rarely spheres. If you do want to sample from a sphere, see the generalization at stats.stackexchange.com/questions/310346. Note that sampling uniformly from a standardized sphere, as described at stats.stackexchange.com/questions/62092/…, and then transforming the sample back to the original units does not give a uniform sample of the ellipsoid. That's why your question still is unclear. $\endgroup$
    – whuber
    Sep 18, 2018 at 21:40
  • 1
    $\begingroup$ I understand that uniformly sampling from the surface of the sphere and then mapping it to ellipsoid won't give uniform samples on the ellipsoid. So I need a method that does sample uniformly from the surface of an ellipsoid. $\endgroup$ Sep 18, 2018 at 22:51
  • 1
    $\begingroup$ Do you want to have the sample uniform on the surface of an ellipsoid, in the sense that each area element dA of the hyper-surface contains the same probability mass? $\endgroup$ Sep 20, 2018 at 11:29

1 Answer 1

5
$\begingroup$

When the different ellipsoid axes are not too much different then it is feasible to use rejection sampling (with large differences you reject a lot making it less feasible)

  • (1) sample on a hyper-sphere
  • (2) squeezing it into a hyper-ellipsoid
  • (3) compute the rate by which the surface area was squeezed
  • (4) reject samples according to that rate.

2D example

example

set.seed(1)
#some matrix to transform n-sphere (in this case 2x2)
m <- matrix(c(1, 0.55, 0.55, 0.55), 2)

# sample multinomial with identity covariance matrix
x <- cbind(rnorm(3000, 0, 1), rnorm(3000, 0, 1))
l1 <- sqrt(x[,1]^2 + x[,2]^2)

# perpendicular vector
per <- cbind(x[,2], -x[,1])

# transform x
x <- x %*% m
# transform perpendicular vector (to see how the area transforms)
per2 <- per %*% m

# get onto unit-"sphere"/ellipsoid
x <- x/l1

# this is how the area contracted
contract <- sqrt(per2[,1]^2 + per2[,2]^2) / sqrt(per[,1]^2 + per[,2]^2)

# then this is how we should choose to reject samples 
p <- contract/max(contract)

# rejecting
choose <- which( rbinom(n=length(p), size=1, p=p) == 1)

#plotting
plot(x[1:length(choose), 1], x[1:length(choose), 2],
     xlim=c(-1.2, 1.2), ylim=c(-1.2, 1.2),
     xlab = expression(x[1]), ylab = expression(x[2]),
     bg=rgb(0, 0, 0, 0.01), cex=0.6, pch=21, col=rgb(0, 0, 0, 0.01))
title("squeezed uniform circle \n ")

#plotting
plot(x[choose,1], x[choose,2],
     xlim=c(-1.2, 1.2), ylim=c(-1.2, 1.2),
     xlab = expression(x[1]), ylab = expression(x[2]),
     bg=rgb(0, 0, 0, 0.01), cex=0.6, pch=21, col=rgb(0, 0, 0, 0.01))
title("squeezed uniform circle \n  with rejection sampling")
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.