Thank you, I changed that and it's much more efficient.
Sigalit.
On 11/13/07, Phil Spector <[EMAIL PROTECTED]> 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
> [EMAIL PROTECTED]
>
>
[[alternative HTML version deleted]]
______________________________________________
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.
