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
[EMAIL PROTECTED] tel:+64(9)373-7599x86730 fax:+64(9)373-7090
www.health.auckland.ac.nz/pharmacology/staff/nholford
S-plus multivariate normal distribution
Nick,
Here is a snippet of code that should get you started.
#
# Set seed and define number of draws desired
#
set.seed(111)
simnum <- 5
#
# Define final parameter estimates from NONMEM output
#
theta <- c(402,-1.11,-11.2,-0.00842, 0.118)
#
# Define variance-covariance matrix of fixed effects estimates
#
varcov <- matrix(NA,length(theta),length(theta))
varcov[,1] <- c(1.34, -1.11E-01, -1.29, -1.83E-04, -1.62E-03)
varcov[,2] <- c(-1.11E-01, 2.50E-01, 2.27E-02, 3.35E-05, 6.37E-04)
varcov[,3] <- c(-1.29, 2.27E-02, 2.10, 8.27E-05, 1.27E-03)
varcov[,4] <- c(-1.83E-04, 3.35E-05, 8.27E-05, 8.69E-07, -4.81E-06)
varcov[,5] <- c(-1.62E-03, 6.37E-04, 1.27E-03, -4.81E-06, 4.00E-04)
varcov <- round(varcov,8)
#
# Generate random draws from MVN for theta values
#
sim.theta <- matrix(NA, nrow = simnum, ncol = length(theta))
sim.theta[,1:length(theta)] <- rmvnorm(simnum,mean = theta, cov =
varcov)
Please let me know if you have other questions.
Steve
Quoted reply history
-----Original Message-----
From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED]
On Behalf Of Nick Holford
Sent: Thursday, September 13, 2007 6:31 AM
To: nmusers
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
[EMAIL PROTECTED] tel:+64(9)373-7599x86730 fax:+64(9)373-7090
www.health.auckland.ac.nz/pharmacology/staff/nholford
Hi Nick,
In R, use mvrnorm (load package MASS first)
In S+, use rmvnorm
Leonid
--------------------------------------
Leonid Gibiansky
President
QuantPharm LLC
www.quantpharm.com
Nick Holford wrote:
> 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
> [EMAIL PROTECTED] tel:+64(9)373-7599x86730 fax:+64(9)373-7090
> www.health.auckland.ac.nz/pharmacology/staff/nholford
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" <[EMAIL PROTECTED]>
Sent by: [EMAIL PROTECTED]
13-Sep-2007 06:30
To
"nmusers" <[email protected]>
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
[EMAIL PROTECTED] tel:+64(9)373-7599x86730 fax:+64(9)373-7090
www.health.auckland.ac.nz/pharmacology/staff/nholford
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
Quoted reply history
-----Original Message-----
From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On
Behalf Of GIRARD PASCAL
Sent: Thursday, September 13, 2007 7:8 AM
To: [EMAIL PROTECTED]; 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
[EMAIL PROTECTED]
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 : [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] De
la part de [EMAIL PROTECTED]
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" <[EMAIL PROTECTED]>
Sent by: [EMAIL PROTECTED]
13-Sep-2007 06:30
To
"nmusers" <[email protected]>
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
[EMAIL PROTECTED] tel:+64(9)373-7599x86730 fax:+64(9)373-7090
www.health.auckland.ac.nz/pharmacology/staff/nholford
Hi,
Thanks very much to all of you who took the time to respond. It was very nice
to see the wide variety of approaches and for the helpful tips for a naive
S-plus user like me.
Nick
Quoted reply history
> -----Original Message-----
> From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED]
> On Behalf Of Nick Holford
> Sent: Thursday, September 13, 2007 6:31 AM
> To: nmusers
> 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
> [EMAIL PROTECTED] tel:+64(9)373-7599x86730 fax:+64(9)373-7090
> www.health.auckland.ac.nz/pharmacology/staff/nholford
--
Nick Holford, Dept Pharmacology & Clinical Pharmacology
University of Auckland, 85 Park Rd, Private Bag 92019, Auckland, New Zealand
[EMAIL PROTECTED] tel:+64(9)373-7599x86730 fax:+64(9)373-7090
www.health.auckland.ac.nz/pharmacology/staff/nholford