I'm trying to generate multivariate data shaped like Saturn (long story).
More formally;
- *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)$).
- 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