>>>>> Martin Krautschke <m.krautsc...@gmx.at>
>>>>> on Mon, 27 Dec 2010 22:47:26 +0100 writes:
> I am looking for an implementation of random rotation
> matrix generation in R to do a rotation test: I want to
> use the matrices to create random multivariate normal
> matrices with common covariance structure and mean based
> on an observed data matrix.
> The rRotationMatrix-function in the mixAK-package is an
> option, but as far as I can tell I need to draw rotation
> matrices with determinant -1 as well. Roast and Romer in
> the limma-bioconductor package appear to have implemented
> something similar, which seems not to be general enough
> for my purposes, however. Inspired by the code in the
> ffmanova-rotationtest function I thought of the following,
> but it appears to me that there only the covariance, not
> the mean, is preserved:
> #####
> # a given Y has independent, multivariate normal rows
> library(mvtnorm)
> Y <- rmvnorm(4,mean=1:10,sigma=diag(1:10)+3)
> # Generation of a set of random matrices Z
> for (i in 1:10) {
> # R is random matrix of independent standard-normal entries
> R <- matrix(rnorm(16),ncol=4)
> R <- qr.Q(qr(R, LAPACK = TRUE))
> # Z shall be a random matrix with the same mean and covariance structure
as Y
> Z <- crossprod(R,Y)
> }
> #####
> A suggestion for the procedure exists (in Dorum et al.
http://www.bepress.com/sagmb/vol8/iss1/art34/ , end of chapter 2.1), but a hint
to a (fast) implementation would be greatly appreciated.
> Best regards and a happy new year,
> Martin Krautschke
> -----------------------
> Martin Krautschke
> Student at University of Vienna
and this is not a home work problem?
Just in case, I don't give you the complete solution in R,
but in words :
Think geometrically:
Rotation in the above sense only preserves the mean when that is
the zero vector.
Consequently: Your procedure must rather be
1) Y0 <- Y - mY
2) Z0 <- Q' %*% Y0
3) Z <- Z0 + mY
and to make this work with data matrices Y, Z,
the mean vector mY must either be a matrix with constant
columns or the result of as.vector()ing such a matrix.
Regards,
Martin Maechler, ETH Zurich
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
