I want to create a random sample of length n from a normal multivariate distribution

I have a problem in understanding this question, especially this part:

"generate a random sample of length n from a normal multivariate"

This is what I have done using the R package mvtnorm:

my_function <- function(n=1,k){
mean = rep(0,k)
sigma = diag(length(mean))
rmvnorm(n, mean, sigma,method=c("eigen", "svd", "chol"), pre0.9_9994 = FALSE)}


my_function(3)

This way the output of my function is, in this example, a vector of three numbers which in my mind is A sample of length 3. Is this correct? Is there another way to do this?

I am asking this question because I have to do the same without using the library rmvnorm and with a normal bivariate.

I just managed to write the density function:

dbivnorm <- function(x,y,mux=0,muy=0,sigmax=1,sigmay=1,rho=0){
(2*pi)^(-1) * ((1-rho^2)*sigmax^2*sigmay^2)^(-.5) *
exp( -((x-mux)^2/(sigmax^2) -2*rho*((x-mux)/sigmax *(y-muy)/sigmay) +
(y-muy)^2/(sigmay^2))/(2*(1-rho^2)) )
}


But now I am not sure how to proceed.

• This is a sample of length one (1) and dimension three (3). Jan 8, 2019 at 18:51
• Ok, how do I get a sample of length n, with let's say n = 3 using rmvnorm? Jan 8, 2019 at 19:23
• Its manual page explains how. You will figure it out quickly enough by running rmvnorm a few times with different argument values. Since (for some reason) you don't want to use rmvnorm, could you please be specific about what tools or techniques you are willing to employ? If the reason for not using the solution available to you is educational, then please tag your post with self-study.
– whuber
Jan 8, 2019 at 19:32
• @whuber: sorry, I typed my answer while you were typing your comment! Jan 8, 2019 at 19:38

Usually, a good starting point to understand how to run an R function is to read the attached documentation:

which also contains an example as an illustration

where the entry n in rmvnorm(n,mean,sigma) is the number of elements in the random sample, while mean is the mean vector and sigma is the variance-covariance matrix.

• Yes, I read that too but I have a doubt on my output: let's say I run rmvnorm(5,mean = rep(0,5), diag(length(mean))) This gives me a 5x5 matrix, but I just want the one random sample of length n, to me this seems 5 random samples of length 5. What am I not getting? Jan 8, 2019 at 19:42
• If $n=5$ and $k=5$, this is a correct outcome of $5$ realisations of a $5$-dimensional random vector. Jan 8, 2019 at 19:50
• In R, a matrix and an array of columns are identical concepts.
– whuber
Jan 8, 2019 at 19:51

Are you allowed to use a package that provides you with independent normal distributed RVs? If not, the Box-Muller method (https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform) could be used for that.

Given that you have n independent normal distributed random variables, using the cholesky decomposition $$L^\top L = \Sigma$$ (for the wanted covariance matrix $$\Sigma$$) and multiplying $$L$$ by the vector of independent normal RVs gives you a multivariate normal distributed random variable with covariance $$\Sigma$$.

• Mmm, could you be a bit more specific on the first part of my questions? I mean, is what I did right for generating a random sample of length n from a normal multivariate? Because for that part I am allowed to use the package. Jan 8, 2019 at 16:00
• I'm very unsure what you are doing in your first part, especially the relation between 'my_function' and 'using_rmvnorm'. Normally, simply 'rmvnorm(rep(0,k),diag(k))' should work.
– nope
Jan 8, 2019 at 16:05
• Yeah, I overcomplicated that, plus using_rmvnorm should just be my_function. Anyway, I have still another question. Given that I have the density function for the bivariate distribution, can't I just use rand() two times to get two values for x and y and then compute dbivnorm(x,y)? Jan 8, 2019 at 16:12