[R] Tricky vectorization problem

[Dimitris Rizopoulos]

a small improvement is to avoid computing the spectral decomposition  
of `Sigma' (look at code of mvrnorm()) in each iteration, e.g.,

#set up a collector for subject means in A
a <- numeric(Ns)
#set up a collector for subject means in B
b <- numeric(Ns)

eS <- eigen(Sigma, symmetric = TRUE, EISPACK = TRUE)
ev <- eS$values
fact <- eS$vectors %*% diag(sqrt(pmax(ev, 0)), 2)

#start a monte carlo experiment
for (i in 1:mce) {
     #generate correlated ideal means for each subject in each condition
     X <- rnorm(2 * Ns)
     dim(X) <- c(2, Ns)
     sub.means <- t(fact %*% X)

     #for each subject
     for (s in 1:Ns) {
         #generate some data for A and grab the mean
         a[s] <- mean(rnorm(a.No[s], sub.means[s, 1], s.a.w[s]))
         #generate some data for B and grab the mean
         b[s] <- mean(rnorm(b.No[s], sub.means[s, 2], s.b.w[s]))
     }
     #store the observed correlation between A and B
     sim.r[i] <- cor(a, b)
}


I hope it helps.

Best,
Dimitris

----
Dimitris Rizopoulos
Ph.D. Student
Biostatistical Centre
School of Public Health
Catholic University of Leuven

Address: Kapucijnenvoer 35, Leuven, Belgium
Tel: +32/(0)16/336899
Fax: +32/(0)16/337015
Web: http://med.kuleuven.be/biostat/
      http://www.student.kuleuven.be/~m0390867/dimitris.htm


Quoting Mike Lawrence <Mike.Lawrence at dal.ca>:

> Seems there were some line break issues when pasting the code, trying
> again with a different commenting style:
>
> #start a timer
> start = proc.time()[1]
>
> #set the true correllation
> rho = .5
>
> #set the number of Subjects
> Ns = 100
>
> #for each subject, set a number of observations in A
> a.No = 1:100
> #for each subject, set a number of observations in B
> b.No = 1:100
>
> #set the between Ss variance in condition A
> v.a = 1
> #set the between Ss variance in condition B
> v.b = 2
>
> #for each subject, set a standard deviation in A
> s.a.w = 1:100
> #for each subject, set a standard deviation in B
> s.b.w = 1:100
>
> #set the number of monte carlo experiments
> mce = 1e3
>
> #set up a collector for the simulated correlations
> sim.r = vector(length=mce)
>
> #define a covariance matrix for use in generating the correlated data
> Sigma=matrix(c(v.a,sqrt(v.a*v.b)*rho,sqrt(v.a*v.b)*rho,v.b),2,2)
>
> #set up a collector for subject means in A
> a = vector(length=Ns)
> #set up a collector for subject means in B
> b = vector(length=Ns)
>
> #start a monte carlo experiment
> for(i in 1:mce){
>
> 	#generate correlated ideal means for for each subject in each condition
> 	sub.means=mvrnorm(Ns,rep(0,2),Sigma)
>
> 	#for each subject
> 	for(s in 1:Ns){
>
> 		#generate some data for A and grab the mean
> 		a[s] = mean(rnorm(a.No[s],sub.means[s,1],s.a.w[s]))
> 		#generate some data for B and grab the mean
> 		b[s] = mean(rnorm(b.No[s],sub.means[s,2],s.b.w[s]))
>
> 	}
>
> 	#store the observed correlation between A and B
> 	sim.r[i] = cor(a,b)
>
> }
>
> #show the total time this took
> print(proc.time()[1]-start)
>
> ______________________________________________
> 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.
>
>



Disclaimer: http://www.kuleuven.be/cwis/email_disclaimer.htm