[R] Re : Elegant Code
Daniel Nordlund
djnordlund at frontier.com
Fri Mar 16 10:04:22 CET 2012
> -----Original Message-----
> From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org]
> On Behalf Of Pascal Oettli
> Sent: Friday, March 16, 2012 12:36 AM
> To: Raphael Fraser
> Cc: r-help at 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.fraser at gmail.com>
> À : r-help at 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
More information about the R-help
mailing list