I am trying to simulate draws from a multivariate normal using mvrnorm, and
am getting the following error message:
Error in mu + eS$vectors %*% diag(sqrt(pmax(ev, 0)), p) %*% t(X) :
non-conformable arrays
I do not understand why I am getting this message, since the vector of means
I am giving to the function is 13 by 1 and the variance matrix I am giving
to the function is 13 by 13...these matrices are therefore conformable,
right?
Is it possible that this is actually a problem with the
positive-definiteness of the variance matrix? The matrix I am attempting to
use has eigenvalues that are close to zero (but positive).
Here's a transcript of the commands I've tried:
> 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
> V
[,1] [,2] [,3] [,4] [,5]
[1,] 2.632324e+05 -1.286920e-01 -9.461830e-02 6.27173e-02 -2.625664e+00
[2,] -1.286920e-01 1.164268e+00 -7.496500e-03 4.95040e-03 -1.557751e-01
[3,] -9.461830e-02 -7.496500e-03 5.761000e-03 -6.66370e-03 6.366440e-02
[4,] 6.271730e-02 4.950400e-03 -6.663700e-03 8.23150e-03 -6.252150e-02
[5,] -2.625664e+00 -1.557751e-01 6.366440e-02 -6.25215e-02 1.753953e+02
[6,] 1.017806e+01 8.797900e-02 -1.028916e+01 9.75542e+00 -8.614497e+03
[7,] -1.286920e-01 1.164268e+00 -7.496500e-03 4.95040e-03 -1.557751e-01
[8,] 3.930000e-07 2.060000e-09 1.150000e-10 -1.39000e-10 5.550000e-09
[9,] -1.950000e-07 3.680000e-09 -2.050000e-11 8.86000e-12 -6.110000e-09
[10,] 4.530000e-07 4.610000e-10 9.250000e-12 -1.50000e-11 1.310000e-08
[11,] 1.170000e-06 -6.710000e-09 6.740000e-11 -1.68000e-10 -3.400000e-09
[12,] -9.600000e-07 3.970000e-09 -9.520000e-11 2.12000e-10 -3.430000e-10
[13,] -4.970000e-07 -2.070000e-10 -3.100000e-11 5.69000e-11 -2.270000e-09
[,6] [,7] [,8] [,9] [,10]
[1,] 1.017806e+01 -1.286920e-01 3.93000e-07 -1.95000e-07 4.53000e-07
[2,] 8.797900e-02 1.164268e+00 2.06000e-09 3.68000e-09 4.61000e-10
[3,] -1.028916e+01 -7.496500e-03 1.15000e-10 -2.05000e-11 9.25000e-12
[4,] 9.755420e+00 4.950400e-03 -1.39000e-10 8.86000e-12 -1.50000e-11
[5,] -8.614497e+03 -1.557751e-01 5.55000e-09 -6.11000e-09 1.31000e-08
[6,] 4.980406e+05 8.797900e-02 -4.67000e-07 3.50000e-07 -7.68000e-07
[7,] 8.797900e-02 1.164268e+00 2.06000e-09 3.68000e-09 4.61000e-10
[8,] -4.670000e-07 2.060000e-09 2.15049e-02 9.00950e-03 9.58560e-03
[9,] 3.500000e-07 3.680000e-09 9.00950e-03 1.58849e-02 1.65230e-03
[10,] -7.680000e-07 4.610000e-10 9.58560e-03 1.65230e-03 1.53024e-02
[11,] 5.610000e-07 -6.710000e-09 4.26370e-03 -9.08500e-04 6.03300e-03
[12,] -2.720000e-07 3.970000e-09 -3.98220e-03 4.55100e-04 -5.55290e-03
[13,] 9.530000e-08 -2.070000e-10 -1.38413e-02 -9.42850e-03 -1.12002e-02
[,11] [,12] [,13]
[1,] 1.17000e-06 -9.60000e-07 -4.97000e-07
[2,] -6.71000e-09 3.97000e-09 -2.07000e-10
[3,] 6.74000e-11 -9.52000e-11 -3.10000e-11
[4,] -1.68000e-10 2.12000e-10 5.69000e-11
[5,] -3.40000e-09 -3.43000e-10 -2.27000e-09
[6,] 5.61000e-07 -2.72000e-07 9.53000e-08
[7,] -6.71000e-09 3.97000e-09 -2.07000e-10
[8,] 4.26370e-03 -3.98220e-03 -1.38413e-02
[9,] -9.08500e-04 4.55100e-04 -9.42850e-03
[10,] 6.03300e-03 -5.55290e-03 -1.12002e-02
[11,] 9.29045e-02 -8.00488e-02 -2.21646e-02
[12,] -8.00488e-02 7.19209e-02 1.82813e-02
[13,] -2.21646e-02 1.82813e-02 1.76818e-02
> B
[,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,] -22.02731 0.4982429 2.651566 -2.356899 31.65277 -3302.719 0.4982429
[,8] [,9] [,10] [,11] [,12] [,13]
[1,] 0.3675237 -0.2562356 0.4549983 4.901309 -4.464062 -0.0836918
> t(B)
[,1]
[1,] -22.0273100
[2,] 0.4982429
[3,] 2.6515660
[4,] -2.3568990
[5,] 31.6527700
[6,] -3302.7190000
[7,] 0.4982429
[8,] 0.3675237
[9,] -0.2562356
[10,] 0.4549983
[11,] 4.9013090
[12,] -4.4640620
[13,] -0.0836918
Thanks very much,
Stu Jordan
Princeton University
[[alternative HTML version deleted]]
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 arraysYou 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. I.e. the point is that> x <- matrix(1:2,1,2) > M <- matrix(1:4,2) > x+MError in x + M : non-conformable arrays> t(x)+MError 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
On 02-Feb-04 Stuart V Jordan wrote:> I am trying to simulate draws from a multivariate normal using mvrnorm, > and > am getting the following error message: > > Error in mu + eS$vectors %*% diag(sqrt(pmax(ev, 0)), p) %*% t(X) : > non-conformable arrays[...] Hmmm ... using the same B and V as giben by Stuart Jordan, I get:> mvrnorm(1,B,V)[1] -179.8332342 0.9632282 2.5687489 -2.2337125 47.2717626 [6] -3745.3844310 0.9632282 0.2839965 -0.1500585 0.5804460 [11] 5.1919420 -4.7381677 -0.2382119> mvrnorm(2,B,V)Error in mu + eS$vectors %*% diag(sqrt(pmax(ev, 0)), p) %*% t(X) : non-conformable arrays ????? Ted. PS By the way, [near] singularity seems to have nothing to do with it here: If I brutally make V non-singular: diag(V) <- diag(V)+10000000 it goes the same way. -------------------------------------------------------------------- E-Mail: (Ted Harding) <Ted.Harding at nessie.mcc.ac.uk> Fax-to-email: +44 (0)870 167 1972 Date: 02-Feb-04 Time: 22:17:13 ------------------------------ XFMail ------------------------------
As Peter Dalgaard has already pointed out, you both need to read the help page. If you supply a matrix where the help page asks for a vector, you are likely to get troubles. I will coerce it for the next version, but I would like to point out that in about 12 years of providing mvrnorm (often as here without credit), this it the first time anyone has made this error in public. Anticipating what people will try is not easy. On Mon, 2 Feb 2004 Ted.Harding at nessie.mcc.ac.uk wrote:> On 02-Feb-04 Stuart V Jordan wrote: > > I am trying to simulate draws from a multivariate normal using mvrnorm, > > and > > am getting the following error message: > > > > Error in mu + eS$vectors %*% diag(sqrt(pmax(ev, 0)), p) %*% t(X) : > > non-conformable arrays > [...] > > Hmmm ... using the same B and V as giben by Stuart Jordan, I get: > > > mvrnorm(1,B,V) > [1] -179.8332342 0.9632282 2.5687489 -2.2337125 47.2717626 > [6] -3745.3844310 0.9632282 0.2839965 -0.1500585 0.5804460 > [11] 5.1919420 -4.7381677 -0.2382119 > > mvrnorm(2,B,V) > Error in mu + eS$vectors %*% diag(sqrt(pmax(ev, 0)), p) %*% t(X) : > non-conformable arrays > > ????? > > Ted. > > PS By the way, [near] singularity seems to have nothing to do > with it here: > If I brutally make V non-singular: > > diag(V) <- diag(V)+10000000 > > it goes the same way. > > > -------------------------------------------------------------------- > E-Mail: (Ted Harding) <Ted.Harding at nessie.mcc.ac.uk> > Fax-to-email: +44 (0)870 167 1972 > Date: 02-Feb-04 Time: 22:17:13 > ------------------------------ XFMail ------------------------------ > > ______________________________________________ > 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 > >-- Brian D. Ripley, ripley at stats.ox.ac.uk Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UK Fax: +44 1865 272595