RE: [NMusers] S-plus multivariate normal distribution

From: Sutjandra, Liviawati <lsutjand_at_amgen.com>
Date: Thu, 13 Sep 2007 15:10:22 -0700


Hi,

Thank you for sharing the code.
It seems that coding with x[row(x)>col(x)] <- x[row(x)<col(x)] does not =
work
properly for matrix bigger than 3x3. Using transpose should work.

x <- matrix(nrow=5, ncol=5)
x[row(x)<=col(x)] <-c(1.00,
                     0.81,1.00,
                     0.32,0.21,1.00,
                     0.11,-0.9,0.03,1.00,
                     -0.05,0.01,0.2,0.51,1.00)
x[row(x)>=col(x)] <- t(x)[row(x)>=col(x)]


Regards,
Liviawati Sutjandra


-----Original Message-----
From: owner-nmusers_at_globomaxnm.com =
[mailto:owner-nmusers_at_globomaxnm.com] On
Behalf Of GIRARD PASCAL
Sent: Thursday, September 13, 2007 7:8 AM
To: jeffrey.a.wald_at_gsk.com; Nick Holford; nmusers
Subject: RE : [NMusers] S-plus multivariate normal distribution

Hi Jeff,

Splus does not like loops because it badly handles RAM, especially if =
you
perform multiple replications. I would rather code the matrix =
initialisation
with:

x <- matrix(nrow=3, ncol=3)
x[row(x)<=col(x)] <-c(1.00,
                     0.81,1.00,
                     0.32,0.21,1.00)
x[row(x)>col(x)] <- x[row(x)<col(x)]

Cheers,

Pascal

Pascal Girard, PhD
EA 3738, CTO
Fac Medecine Lyon-Sud, BP12
69921 OULLINS Cedex, France
pascal.girard_at_adm.univ-lyon1.fr
Tel  +33 (0)4 26 23 59 54 / Fax +33 (0)4 26 23 59 76
 
Master Recherche Lyon 1 Santé et Populations, Spécialité =
PhIT 
http://master-sante-pop.univ-lyon1.fr/
-----Message d'origine-----
De : owner-nmusers_at_globomaxnm.com =
[mailto:owner-nmusers_at_globomaxnm.com] De
la part de jeffrey.a.wald_at_gsk.com
Envoyé : jeudi 13 septembre 2007 15:40
À : Nick Holford; nmusers
Objet : Re: [NMusers] S-plus multivariate normal distribution


Nick - Here is an adapted snippet of some code I have used recently.
 Implementation of multivariate normals is pretty easy in S+, but it =
can be
an annoying pain in the tuckus  to take the lower half of the =
covariance
matrix from NONMEM and create the symmetrical matrix that S is =
expecting
(especially when you have a big matrix).  There is some code, =
admittedly
awkward, that will take NONMEM's half matrix (prettied up with some =
commas)
and create a symmetrical matrix.  In this case I am using a =
correlation
matrix, but you'll get the idea.  

Regards, Jeff

nsim<-100

### A correlation matrix
a<-c(1.00,
         0.81,1.00,
         0.32,0.21,1.00)

a<-matrix(a, ncol=1, byrow=T)

ndim <- 3
sel<-0
# Will make symetric matrix from the lower triangle called b
b<-matrix(0,ncol=ndim,nrow=ndim)
for (i in seq(ndim)) {
  for (j in seq(i)) {
    sel <- sel+1
    b[i,j]<-a[sel]}}
for (i in seq(ndim-1)) {
  for (j in seq(ndim-i)) {
    sel <- j+i
    b[i,sel]<-b[sel,i]}}

stuff<-rmvnorm(nsim, mean=c(0,0,0), cov=b, sd=c(0.15, 0.30,0.5))

# make a plot
splom(~stuff)

#-##-##-##-##-##-##-##-##-##-##-##-##-##-##-##-###-#




"Nick Holford" <n.holford_at_auckland.ac.nz>
Sent by: owner-nmusers_at_globomaxnm.com
13-Sep-2007 06:30
       
To
"nmusers" <nmusers_at_globomaxnm.com>
cc

Subject
[NMusers] S-plus multivariate normal distribution


Hi,

Can someone give me some clues on how to sample from a multivariate =
normal
distribution using S-plus?

Some working code examples would be very helpful.

Thanks,

Nick

--
Nick Holford, Dept Pharmacology & Clinical Pharmacology
University of Auckland, 85 Park Rd, Private Bag 92019, Auckland, New =
Zealand
n.holford_at_auckland.ac.nz tel:+64(9)373-7599x86730 fax:+64(9)373-7090
www.health.auckland.ac.nz/pharmacology/staff/nholford

Received on Thu Sep 13 2007 - 18:10:22 EDT

This archive was generated by hypermail 2.2.0 : Tue Nov 06 2007 - 15:07:18 EST