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

This may sound like a stupid question, but 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). – Xi'an Jan 8 at 18:51
• Ok, how do I get a sample of length n, with let's say n = 3 using rmvnorm? – qcc101 Jan 8 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 at 19:32
• @whuber: sorry, I typed my answer while you were typing your comment! – Xi'an Jan 8 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? – qcc101 Jan 8 at 19:42
• If $n=5$ and $k=5$, this is a correct outcome of $5$ realisations of a $5$-dimensional random vector. – Xi'an Jan 8 at 19:50
• In R, a matrix and an array of columns are identical concepts. – whuber Jan 8 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. – qcc101 Jan 8 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 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)? – qcc101 Jan 8 at 16:12