[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