[R] mvrnorm problem

[Torsten Hothorn]

> Stuart V Jordan <sjordan at princeton.edu> writes:
>
> > > mvrnorm(n = 1000,B,V)
> > Error in mu + eS$vectors %*% diag(sqrt(pmax(ev, 0)), p) %*% t(X) :
> >         non-conformable arrays
> > > mvrnorm(n = 1000,t(B),V)
> > Error in mu + eS$vectors %*% diag(sqrt(pmax(ev, 0)), p) %*% t(X) :
> >         non-conformable arrays
>
> You might, for at least two good reasons, have said that this is from
> library(MASS). The point is that
>
> > mvrnorm(n=10,matrix(c(1,1),1,2),diag(2))
> Error in mu + eS$vectors %*% diag(sqrt(pmax(ev, 0)), p) %*% t(X) :
>         non-conformable arrays
> > mvrnorm(n=10,matrix(c(1,1),2,1),diag(2))
> Error in mu + eS$vectors %*% diag(sqrt(pmax(ev, 0)), p) %*% t(X) :
>         non-conformable arrays
> > mvrnorm(n=10,c(1,1),diag(2))
>             [,1]       [,2]
>  [1,]  0.5005327  1.1919216
>  [2,]  2.8273925  2.7004788
>  [3,]  2.6493970  1.1304274
> ....
>
> and the docs quite clearly say that mu wants to be a vector, not a
> matrix.
>
> Curiously enough, this works with rmvnorm from the mvtnorm package by
> Genz, Bretz, and Hothorn, the difference being that this version adds
> in the means with a sweep() operation, whereas mvrnorm just adds mu
> (to the transpose of the ultimate result) and relies on recycling
> rules.

Credits to Fritz: {rd}mvnorm moved from `e1071' to `mvtnorm' for
obvious reasons some time ago.

Best,

Torsten

> I.e. the point is that
>
> > x <- matrix(1:2,1,2)
> > M <- matrix(1:4,2)
> > x+M
> Error in x + M : non-conformable arrays
> > t(x)+M
> Error in t(x) + M : non-conformable arrays
> > c(x)+M
>      [,1] [,2]
> [1,]    2    4
> [2,]    4    6
> > sweep(M,1,x,"+")
>      [,1] [,2]
> [1,]    2    4
> [2,]    4    6
>
>
> --
>    O__  ---- Peter Dalgaard             Blegdamsvej 3
>   c/ /'_ --- Dept. of Biostatistics     2200 Cph. N
>  (*) \(*) -- University of Copenhagen   Denmark      Ph: (+45) 35327918
> ~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk)             FAX: (+45) 35327907
>
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://www.stat.math.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
>
>