On 07-12-2012, at 18:01, arun wrote: > > > HI, > > I can see ?runif(), ?rnorm() without a set seed.
Good eyesight! Correct. Oops. Inserting set.seed(413) at the start I now get this # > benchmark(eta1 <- f1(X, B, x, sem1), eta2 <- f2(X, B, x, sem1), eta3 <- f3(X, B, x, sem1), # + eta4 <- f4(X, B, x, sem1), eta5 <- f5(X, B, x, sem1), eta6 <- f6(X, B, x, sem1), # + replications=10, columns=c("test","elapsed","relative")) # test elapsed relative # 1 eta1 <- f1(X, B, x, sem1) 2.248 1.392 # 2 eta2 <- f2(X, B, x, sem1) 1.738 1.076 # 3 eta3 <- f3(X, B, x, sem1) 2.024 1.253 # 4 eta4 <- f4(X, B, x, sem1) 2.072 1.283 # 5 eta5 <- f5(X, B, x, sem1) 2.038 1.262 # 6 eta6 <- f6(X, B, x, sem1) 1.615 1.000 Relative conclusions stay the same. The compiled versions f2, f6 are the quickest. Berend > A.K. > > ----- Original Message ----- > From: David L Carlson <dcarl...@tamu.edu> > To: 'arun' <smartpink...@yahoo.com>; "'Doran, Harold'" <hdo...@air.org> > Cc: 'R help' <r-help@r-project.org>; 'David Winsemius' > <dwinsem...@comcast.net> > Sent: Friday, December 7, 2012 11:54 AM > Subject: RE: [R] Vectorizing integrate() > > Must be a cut and paste issue. All three agree on the results but they are > different from those in arun's message: > >> B <- c(0,1) >> sem1 = runif(10, 1, 2) >> x <- rnorm(10) >> X <- cbind(1, x) >> eta <- numeric(10) >> for(j in 1:nrow(X)){ > + fun <- function(u) 1/ (1 + exp(- (B[1] + B[2] * (x[j] + u)))) * dnorm(u, > 0, sem1[j]) > + eta[j] <- integrate(fun, -Inf, Inf)$value > + } >> eta > [1] 0.6586459 0.5565611 0.5226462 0.6824012 0.6908079 0.7578524 0.4715976 > [8] 0.2165210 0.2378657 0.3492133 >> fun <- function(u, m, s) 1/ (1 + exp(- (B[1] + B[2] * > + (m + u)))) * dnorm(u, 0, s) >> eta2 <- sapply(1:nrow(X), function(i) integrate(fun, > + -Inf, Inf, m=x[i], s=sem1[i])$value) >> eta2 > [1] 0.6586459 0.5565611 0.5226462 0.6824012 0.6908079 0.7578524 0.4715976 > [8] 0.2165210 0.2378657 0.3492133 >> identical(eta, eta2) > [1] TRUE >> res<-mapply(function(i) > integrate(fun,-Inf,Inf,m=x[i],s=sem1[i])$value,1:nrow(X)) >> res > [1] 0.6586459 0.5565611 0.5226462 0.6824012 0.6908079 0.7578524 0.4715976 > [8] 0.2165210 0.2378657 0.3492133 >> identical(eta, res) > [1] TRUE > > ------- > David C > >> -----Original Message----- >> From: arun [mailto:smartpink...@yahoo.com] >> Sent: Friday, December 07, 2012 10:36 AM >> To: Doran, Harold >> Cc: R help; David L Carlson; David Winsemius >> Subject: Re: [R] Vectorizing integrate() >> >> >> >> Hi, >> Using David's function: >> fun <- function(u, m, s) 1/ (1 + exp(- (B[1] + B[2] * >> (m + u)))) * dnorm(u, 0, s) >> >> res<-mapply(function(i) integrate(fun,- >> Inf,Inf,m=x[i],s=sem1[i])$value,1:nrow(X)) >> res >> # [1] 0.5212356 0.6214989 0.5306124 0.5789414 0.3429795 0.6972879 >> 0.5952949 >> #[8] 0.7531899 0.4740685 0.7576587 >> identical(res,eta) >> #[1] TRUE >> A.K. >> >> ----- Original Message ----- >> From: "Doran, Harold" <hdo...@air.org> >> To: David Winsemius <dwinsem...@comcast.net> >> Cc: "r-help@r-project.org" <r-help@r-project.org> >> Sent: Friday, December 7, 2012 10:14 AM >> Subject: Re: [R] Vectorizing integrate() >> >> David et al >> >> Thanks, I should have made the post more complete. I routinely use >> apply functions, but often avoid mapply() as I find it so non- >> intuitive. In this instance, I think the situation demands I change >> that position, though. >> >> Reproducible code for the current implementation of the function is >> >> B <- c(0,1) >> sem1 = runif(10, 1, 2) >> x <- rnorm(10) >> X <- cbind(1, x) >> eta <- numeric(10) >> >> for(j in 1:nrow(X)){ >> fun <- function(u) 1/ (1 + exp(- (B[1] + B[2] * (x[j] + u)))) * >> dnorm(u, 0, sem1[j]) >> eta[j] <- integrate(fun, -Inf, Inf)$value >> } >> >> I can't get my head around how mapply() would work here. It accepts as >> its first argument a function. But, in my case I have two functions: >> the user defined integrand, fun(), an then of course calling the R >> function integrate(). >> >> I was thinking maybe along these lines, but this is obviously wrong. >> >> mapply(integrate(function(u) 1/ (1 + exp(- (B[1] + B[2] * (x + u)))) * >> dnorm(u, 0, sem1), -Inf, Inf)$value, MoreArgs = list(B, x, sem1)) >> >> >> >>> -----Original Message----- >>> From: David Winsemius [mailto:dwinsem...@comcast.net] >>> Sent: Thursday, December 06, 2012 1:59 PM >>> To: Doran, Harold >>> Cc: r-help@r-project.org >>> Subject: Re: [R] Vectorizing integrate() >>> >>> >>> On Dec 6, 2012, at 10:10 AM, Doran, Harold wrote: >>> >>>> I have written a program to solve a particular logistic regression >> problem >>> using IRLS. In one step, I need to integrate something out of the >> linear >>> predictor. The way I'm doing it now is within a loop and it is as you >> would >>> expect slow to process, especially inside an iterative algorithm. >>>> >>>> I'm hoping there is a way this can be vectorized, but I have not >> found >>>> it so far. The portion of code I'd like to vectorize is this >>>> >>>> for(j in 1:nrow(X)){ >>>> fun <- function(u) 1/ (1 + exp(- (B[1] + B[2] * (x[j] + u)))) * >> dnorm(u, 0, >>> sd[j]) >>>> eta[j] <- integrate(fun, -Inf, Inf)$value } >>>> >>> >>> The Vectorize function is just a wrapper to mapply. If yoou are able >> to get >>> that code to execute properly for your un-posted test cases, then why >> not >>> use mapply? >>> >>> >>>> Here X is an n x p model matrix for the fixed effects, B is a >> vector with the >>> estimates of the fixed effects at iteration t, x is a predictor >> variable in the jth >>> row of X, and sd is a variable corresponding to x[j]. >>>> >>>> Is there a way this can be done without looping over the rows of X? >>>> >>>> Thanks, >>>> Harold >>>> >>>> [[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. >>> >>> David Winsemius, MD >>> Alameda, CA, 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. > > ______________________________________________ > 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. ______________________________________________ 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.