5
$\begingroup$

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 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
$\endgroup$
  • $\begingroup$ What do you mean exactly by "where every out of diagonal entry of R is the same number in (0,1)"? Can you give an example, such as for p=3? $\endgroup$ – petrichor Jun 16 '11 at 11:17
  • $\begingroup$ yes:see update3 please. $\endgroup$ – user603 Jun 16 '11 at 13:51
4
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ almost. But i'm struggling with a small implementation issue, see edit to the question. $\endgroup$ – user603 Jun 16 '11 at 12:34
  • $\begingroup$ Okay, but then don't two Gaussians with diagonal covariance do the trick (making sure that the ordering of the eigenvalues matches up, and then scaling the "rings")? $\endgroup$ – Bob Durrant Jun 16 '11 at 12:47
  • $\begingroup$ humm -- that is precisely my problem, the scaling. See update2 $\endgroup$ – user603 Jun 16 '11 at 13:33

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.