How to interpret multiple calls of rnorm() function in R?

Question

I'm using rnorm() to generate data for a trainng set. The target is to generate a matrix X with n rows and p columns to represent n sets of p features, and I did so in two ways:

First,

X = matrix(rnorm(n*p), n, p)

Second,

X = matrix(ncol = p, nrow = n)
for (i in 1:p) {
  X[, i] = rnorm(n)
}

I was expecting both outputs to be i.i.d. sample data from a normal distribution with mean=0 and sd=1. But I got confused to see that the model coefficients fitted on these two data sets are very different.

What am I missing here? Thanks.

Answer 1

This question is really about R syntax and not about the normal distribution or sampling. So I guess it's 'off topic' here, and may be closed.

I'm not quite sure what your difficulty is. If you set the same seed before each program, then you should fill an $n\times p$ matrix by columns with exactly the same normal variates.

n=5; p=2
set.seed(713)            # set seed
X = matrix(rnorm(n*p), n, p); X
           [,1]        [,2]
[1,] -1.5925966  0.37590371
[2,] -0.1458793  0.28096447
[3,] -1.0048927 -0.60433386
[4,] -1.5337787 -0.03059243
[5,]  0.3318419  0.46158609
set.seed(713)           # set same seed again
X = matrix(ncol = p, nrow = n)
for (i in 1:p) {
 X[, i] = rnorm(n) }
X
           [,1]        [,2]
[1,] -1.5925966  0.37590371
[2,] -0.1458793  0.28096447
[3,] -1.0048927 -0.60433386
[4,] -1.5337787 -0.03059243
[5,]  0.3318419  0.46158609

---

But if I've missed the point and you're asking something else, then I suggest you replace randomly generated normal variates by sequences of integers. Then the rules for filling matrices might be clearer to you:

n = 5; p = 2; matrix(1:(n*p),n,p)  # n rows, p col
     [,1] [,2]
[1,]    1    6
[2,]    2    7
[3,]    3    8
[4,]    4    9
[5,]    5   10

n = 5; p = 2
X = matrix(ncol = p, nrow = n)
for (i in 1:p) {
 X[, i] = 1:n
}
X
     [,1] [,2]
[1,]    1    1
[2,]    2    2
[3,]    3    3
[4,]    4    4
[5,]    5    5

When matrix turns a vector into a matrix, the default order is to fill the matrix by columns. (All vectors, unless specifically modified, are considered column vectors, even if they print out as rows.)

y = 1:10;  matrix(y, ncol=2)
     [,1] [,2]
[1,]    1    6
[2,]    2    7
[3,]    3    8
[4,]    4    9
[5,]    5   10

y = 1:10;  matrix(y, ncol=2, byrow=T)
     [,1] [,2]
[1,]    1    2
[2,]    3    4
[3,]    5    6
[4,]    7    8
[5,]    9   10

You can make row vectors.

w = c(1,2,3,4,3,2)   # 'c' for column vector
X = as.matrix(w)
X
     [,1]
[1,]    1
[2,]    2
[3,]    3
[4,]    4
[5,]    3
[6,]    2

t(X)  # transpose is a row vector
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]    1    2    3    4    3    2

t(X) %*% X   # matrix multiplication
     [,1]
[1,]   43
X %*% t(X) 
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]    1    2    3    4    3    2
[2,]    2    4    6    8    6    4
[3,]    3    6    9   12    9    6
[4,]    4    8   12   16   12    8
[5,]    3    6    9   12    9    6
[6,]    2    4    6    8    6    4
X*t(X)
Error in X * t(X) : non-conformable arrays
t(w)
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]    1    2    3    4    3    2
t(w)%*%w
     [,1]
[1,]   43

Answer 2

The question is about R syntax, as BruceET says, but I think that data simulation is an important topic, and that you could (should?) generate your data in another way. You can use the mvtnorm package to generate random multivariate matrices and the scale() function to ensure that each column has mean=0 and variance=1:

> library(mvtnorm)
> n <- 6
> p <- 3
> mean <- rep(0, p)
> sigma <- diag(p)              # identity matrix
> X <- rmvnorm(n, mean=mean, sigma=sigma)
> X
           [,1]       [,2]       [,3]
[1,] -0.4557803 -1.6174866 0.34224496
[2,]  0.7982673 -0.1214584 3.05775038
[3,] -0.2175150  1.1449116 0.07709425
[4,]  0.9070137  0.4484822 1.49516553
[5,] -1.8435512 -2.4068282 1.23142658
[6,] -0.5840670 -1.3760948 0.37929271
> apply(X, MARGIN=2, FUN=mean)
[1] -0.2326054 -0.6547457  1.0971624
> apply(X, MARGIN=2, FUN=sd)
[1] 1.012930 1.360682 1.108703
> X <- scale(X)
> apply(X, MARGIN=2, FUN=mean)
[1]  1.127570e-17 -3.699840e-17  9.483155e-17
> apply(X, MARGIN=2, FUN=sd)
[1] 1 1 1

The columns are correlated when $n$ is small, but the correlation decreases (as expected, because sigma is an identity matrix) as $n$ grows:

> X <- scale(rmvnorm(1000, mean=mean, sigma=sigma))
> cor(X)
            [,1]        [,2]        [,3]
[1,]  1.00000000 -0.05327000 -0.01848098
[2,] -0.05327000  1.00000000 -0.01011558
[3,] -0.01848098 -0.01011558  1.00000000
> X <- scale(rmvnorm(10000, mean=mean, sigma=sigma))
> cor(X)
            [,1]        [,2]        [,3]
[1,] 1.000000000 0.005957725 0.002865598
[2,] 0.005957725 1.000000000 0.008932789
[3,] 0.002865598 0.008932789 1.000000000

You can generate heteroscedastic and corrrelated data by changing sigma.