How to generate a pair of ellipses shaped as Saturn in $\mathbb{R}^p$

Question

I'm trying to generate multivariate data shaped like Saturn (long story).

More formally;

1. *cluster 1 is a rank-p Gaussian with correlation matrix R (where every out of diagonal entry of R is the same number in $(0,1)$).
2. cluster 2 is a bunch of points distributed on a hyper-plan of rank p-1, around the equator of cluster 1, all located beyond a certain mahalanobis distance wrt to cluster1 of say $\zeta$.

Rejection sampling is ok, so long as the rejection rate is small: this configuration has to be generated many times for $n$ and $p$ large reliably.

Update: So i implemented Bob Durrant's solution:

x0<-matrix(rnorm(n*p),n,p)
b0<-matrix(rnorm(n*p),n,p)
b1<-qchisq(0.99,df=p);
b2<-sqrt(c(b1,b1*1.25))
b0<-b0/sqrt(rowSums(b0*b0))*runif(n,b2[1],b2[2])
plot(rbind(x0,x1))
mahalanobis(cbind(0,b0),colMeans(x0),var(x0))/(qchisq(0.99,df=p))

yielding the attached picture, which look indeed like Saturn . with the rings located at least @ qchisq(0.99,df=p) away from the center.

Now however, i want my saturn to be ellipse shaped, i.e. to have correlation structure R. The problem is that is i pre-multiply $\verb+x0+$ and $\verb+cbind(0,x1)+$ by $R^{1/2}$, the rings are no longer on the equator :(

Update2:

For example, things already go awry when i replace the diagonal variance structure above with something else.

For example:

library(MASS)
x0<-mvrnorm(n,rep(0,p),diag(rchisq(p,p),p))
b0<-mvrnorm(n,rep(0,p-1),diag(rchisq(p-1,p-1),p-1))
b1<-qchisq(0.99,df=p);
b2<-sqrt(c(b1,b1*1.25))
b0<-b0/sqrt(rowSums(b0*b0))*runif(n,b2[1],b2[2])

mahalanobis(cbind(0,b0),colMeans(x0),var(x0))/(qchisq(0.99,df=p))

most of the distances are way too large (compare with the output to the same call after using diagonal co-variance matrix above)!

Update3

delta<-0.9
p<-3
R<-matrix(runif(p^2,delta*0.99,delta),p,p) 
#to avoid repeating eigen-values, i jitter a bit.
diag(R)<-1

Answer 1

For large $p$ a spherical rank $p$ Gaussian in $\mathbb{R}^{p}$ looks like the uniform distribution on the surface of the hypersphere $\mathbb{S}^{p-1}$ with radius $\sigma\sqrt{p}$, while a rank $p-1$ spherical Gaussian embedded in $\mathbb{R}^{p}$ will give Saturn's rings (the points will look like the hypersphere $\mathbb{S}^{p-2}$).

So I think you can generate this data by drawing from two spherical Gaussians $\mathcal{N}(0,\sigma_{p}^{2}I_{p})$ and $\mathcal{N}(0,\sigma_{p-1}^{2}I_{p-1})$ having set $\sigma_{p}^{2}$ and $\sigma_{p-1}^{2}$ to get the separation you want.

The concentration in the norms is exponentially fast w.r.t $p$, so you probably won't have to throw away any points if $p$ is large enough.