12
$\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$
  • 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 '18 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 '18 at 14:07
  • 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$ – sachin vernekar Sep 18 '18 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$ – Martijn Weterings Sep 20 '18 at 11:29
  • 1
    $\begingroup$ Why, how and where are you going to apply this uniform sample? Such information may help to come with a best/sufficient strategy. For instance, when the different ellipsoid axes are not to much different then you can use rejection sampling by (1) sampling on a sphere, (2) squeezing it into an ellipsoid, (3) compute the rate by which the surface area was squeezed (4) reject samples according to the inverse of that rate. $\endgroup$ – Martijn Weterings Sep 20 '18 at 11:54
3
$\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, privacy policy and cookie policy

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