Thank you, I changed that and it's much more efficient.
Sigalit.


On 11/13/07, Phil Spector <spector@stat.berkeley.edu> wrote:
>
> You can use the replicate function to do your simulation.  First,
> put the code to do one repetition in a function:
>
> dosim0 = function(n=500){
>   x <- 0
>   y <- 0
>   z <- 0
>   aps <- 0
>   tiss <- 0
>   for (i in 1:n){
>      z[i] <- rbinom(1, 1, .6)
>      x[i] <- rbinom(1, 1, .95)
>      y[i] <- z[i]*x[i]
>      if (y[i]==1) aps[i] <- rnorm(1,mean=13.4, sd=7.09) else aps[i] <-
> rnorm(1,mean=12.67, sd=6.82)
>      if (y[i]==1) tiss[i] <- rnorm(1,mean=20.731,sd=9.751) else  tiss[i]
> <- rnorm(1,mean=18.531,sd=9.499)
>    }
>    v <- data.frame(y, aps, tiss)
>    glm(y~., family=binomial, data=v)$coef
> }
>
>
> Now you can run
>
> results = t(replicate(1000,dosim0()))
>
> to get a 1000x3 matrix with the coefficients from each simulation.
>
> Before you actually do the simulations, you might want to rewrite
> your function so that it's not so inefficient.  Functions in R are
> vectorized, and you should take advantage of that whenever possible.
> Here's a vectorized version of your simulation function:
>
> dosim = function(n=500){
>    z = rbinom(n,1,.6)
>    x = rbinom(n,1,.95)
>    y = z * x
>    aps = ifelse(y == 1,rnorm(n,mean=13.4,sd=7.09),rnorm(n,mean=12.67, sd=
> 6.82))
>    tiss = ifelse(y == 1,rnorm(n,mean=20.731,sd=9.751),rnorm(n,mean=18.531
> ,sd=9.499))
>    glm(y~.,family=binomial,data=data.frame(y,aps,tiss))$coef
> }
>
> Now, you can run
>
> results = t(replicate(1000,dosim()))
>
> to get a similar result.  Notice the difference in time for the
> two functions:
>
> > system.time(results <- t(replicate(1000,dosim0())))
>    user  system elapsed
> 55.123   0.004  55.239
> > system.time(results <- t(replicate(1000,dosim())))
>    user  system elapsed
> 20.297   0.000  20.295
>
>                                        - Phil Spector
>                                         Statistical Computing Facility
>                                         Department of Statistics
>                                         UC Berkeley
>                                         spector@stat.berkeley.edu
>
>

	[[alternative HTML version deleted]]

