[R] help about random generation of a Normal distribution of

(Ted Harding) Ted.Harding at manchester.ac.uk
Sat Jul 5 19:09:29 CEST 2008


There is no need to load the MASS library, since the code for
mvrnorm therein is compact and self-contained:

mvrnorm <- function (n=1, mu, Sigma, tol=1e-06, empirical=FALSE) 
{
  p <- length(mu)
  if(!all(dim(Sigma) == c(p, p))) 
    stop("incompatible arguments")
  eS <- eigen(Sigma, symmetric = TRUE, EISPACK = TRUE)
  ev <- eS$values
  if(!all(ev >= -tol * abs(ev[1]))) 
    stop("'Sigma' is not positive definite")
  X <- matrix(rnorm(p * n), n)
  if(empirical) {
    X <- scale(X, TRUE, FALSE)
    X <- X %*% svd(X, nu = 0)$v
    X <- scale(X, FALSE, TRUE)
  }
  X <- drop(mu) +
       eS$vectors %*% diag(sqrt(pmax(ev, 0)), p) %*% t(X)
  nm <- names(mu)
  if(is.null(nm) && !is.null(dn <- dimnames(Sigma))) 
    nm <- dn[[1]]
  dimnames(X) <- list(nm, NULL)
  if(n == 1) 
    drop(X)
  else t(X)
}


Define that function as above, then proceed along the lines suggested
by Gavin Simpson below.

Ted.


On 05-Jul-08 16:43:46, Gavin Simpson wrote:
> On Sat, 2008-07-05 at 18:21 +0200, Arnau Mir wrote:
>> Hello.
>> 
>> Somebody knows how can I generate a set of n random vectors of a
>> normal distribution of several variables?
>> For example, I want to generate n=100 random vectors of two
>> dimensions for a normal with mean c(0,1) and variance matrix: 
>> matrix(c(2,1,1,3),2,2).
> 
> One is mvrnorm() in the MASS package, part of the VR bundle that comes
> with R.
> 
>> require(MASS)
>> mu <- c(0,1)
>> Sigma <- matrix(c(2,1,1,3),2,2)
>> res <- mvrnorm(100, mu = mu, Sigma = Sigma)
>> head(res)
>            [,1]        [,2]
> [1,]  2.7582876  1.04208798
> [2,]  0.6364184 -0.08043244
> [3,] -1.8897731  0.04051395
> [4,]  2.6699881  0.83163661
> [5,] -1.1942385 -1.17503716
> [6,] -0.4303459 -0.80880649
> 
> HTH
> 
> G

--------------------------------------------------------------------
E-Mail: (Ted Harding) <Ted.Harding at manchester.ac.uk>
Fax-to-email: +44 (0)870 094 0861
Date: 05-Jul-08                                       Time: 18:09:23
------------------------------ XFMail ------------------------------



More information about the R-help mailing list