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

[(Ted Harding)]

Bill, while your point is a fair one, as far as it goes,
I (and no doubt others) often run R on low-powered hardware,
and therefore like to restrict overload. MASS is quite a big
package.

I've no dispute with your statement that MASS is both recommended
and universally available. And it is valuable and good.I was simply
pointing out that, if all you want from MASS is mvrnorm, then you
don't need to load MASS, once you have extracted mvrnorm's definition
from MASS.

BUT [OVERSIGHT IN MY FIRST POST]: As well as the code for mvrnorm,
you must also define the ginv() function (needed by mvrnorm):

ginv <- function (X, tol = sqrt(.Machine$double.eps)) 
{
  if(length(dim(X)) > 2 || !(is.numeric(X) || is.complex(X))) 
    stop("'X' must be a numeric or complex matrix")
  if (!is.matrix(X)) X <- as.matrix(X)
  Xsvd <- svd(X)
  if (is.complex(X)) Xsvd$u <- Conj(Xsvd$u)
  Positive <- Xsvd$d > max(tol * Xsvd$d[1], 0)
  if (all(Positive)) Xsvd$v %*% (1/Xsvd$d * t(Xsvd$u))
    else if (!any(Positive)) array(0, dim(X)[2:1])
      else Xsvd$v[, Positive, drop = FALSE] %*% ((1/Xsvd$d[Positive]) * 
        t(Xsvd$u[, Positive, drop = FALSE]))
}


For what it's worth, I use mvrnorm so often (also ginv) that the code
for mvrnorm (as below) and for ginv (as above) is part of my .Rprofile,
so they are always available without having to load anything else.

Best wishes,
Ted.

On 06-Jul-08 05:36:12, Bill.Venables at csiro.au wrote:
> What's the problem with loading the MASS library???  
> 
> The MASS library is one of the recommended packages, so should be
> universally available.
> 
> Pinching code like this and freezing it inside your personal scripts
> means that if ever we find bug fixes or improvements you miss out on
> them. 
> 
> 
> Bill Venables
> CSIRO Laboratories
> AUSTRALIA
> http://www.cmis.csiro.au/bill.venables/ 
> 
> -----Original Message-----
> From: r-help-bounces at r-project.org
> [mailto:r-help-bounces at r-project.org] On Behalf Of Peng Jiang
> Sent: Sunday, 6 July 2008 1:47 PM
> To: Arnau Mir Torres
> Cc: r-help at r-project.org
> Subject: Re: [R] help about random generation of a Normal distribution
> ofseveral variables
> 
> Hi , Arnau
> 
> Did you ever check your mailbox?  your question was answered last  
> night  Beijing time. :)
> Just read the following .
> -------------------------
> 
> 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 ------------------------------
> 
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
> On 2008-7-6, at ä¸Å_午12:13, Arnau Mir Torres 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).
>>
>> Thanks in advance,
>>
>> Arnau.
>>
>> ______________________________________________
>> R-help at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide
>> http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
> 
>   
>       [[alternative HTML version deleted]]
> 
> 
> 
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

--------------------------------------------------------------------
E-Mail: (Ted Harding) <Ted.Harding at manchester.ac.uk>
Fax-to-email: +44 (0)870 094 0861
Date: 06-Jul-08                                       Time: 11:57:29
------------------------------ XFMail ------------------------------