> -----Original Message----- > From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] > On Behalf Of Pascal Oettli > Sent: Friday, March 16, 2012 12:36 AM > To: Raphael Fraser > Cc: r-help@r-project.org > Subject: [R] Re : Elegant Code > > Hi Raphael, > > Something like that? > > require(pscl)sim <- numeric() > for(i in 1:5){ > eval(parse(text=paste('b',i,' <- rigamma(50,1,1)',sep=''))) > eval(parse(text=paste('theta',i,' <- > rgamma(50,0.5,(1/b',i,'))',sep=''))) > eval(parse(text=paste('sim',i,' <- rpois(50,theta',1,')',sep=''))) > eval(parse(text=paste('sim <- cbind(sim, sim',i,')',sep=''))) > } > x11() > par(mfrow=c(1,5)) > boxplot(c(sim1,sim2,sim3,sim4,sim5)) > > x11() > boxplot(sim) > > > Regards, > Pascal > > > ----- Mail original ----- > De : Raphael Fraser <raphael.fra...@gmail.com> > À : r-help@r-project.org > Cc : > Envoyé le : Vendredi 16 mars 2012 16h09 > Objet : [R] Elegant Code > > Hi, > > Can anyone help to write a more elegant version of my code? I am sure > this can be put into a loop but I am having trouble creating the > objects b1,b2,b3,...,etc. > > b1 <- rigamma(50,1,1) > theta1 <- rgamma(50,0.5,(1/b1)) > sim1 <- rpois(50,theta1) > > b2 <- rigamma(50,1,1) > theta2 <- rgamma(50,0.5,(1/b2)) > sim2 <- rpois(50,theta2) > > b3 <- rigamma(50,1,1) > theta3 <- rgamma(50,0.5,(1/b3)) > sim3 <- rpois(50,theta3) > > b4 <- rigamma(50,1,1) > theta4 <- rgamma(50,0.5,(1/b4)) > sim4 <- rpois(50,theta4) > > b5 <- rigamma(50,1,1) > theta5 <- rgamma(50,0.5,(1/b5)) > sim5 <- rpois(50,theta5) > > > > par(mfrow=c(1,5)) > boxplot(sim1) > boxplot(sim2) > boxplot(sim3) > boxplot(sim4) > boxplot(sim5); > > Thanks, > Raphael >
Or, just get rid of the for loops altogether, something like this trials <-5 N <- 50 b <- rigamma(N*trials, 1, 1) theta <- rgamma(N*trials,0.5,1/b) sim <- rpois(N*trials,theta) dim(sim) <- c(N,trials) boxplot(sim) Hope this is helpful, Dan Daniel Nordlund Bothell, WA USA ______________________________________________ R-help@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.