vectorizing the integrate function

10 messages
Open this post in threaded view
|

vectorizing the integrate function

 Hi all,I am having some difficulties in vectorizing the integrate function. Let me explain with an example. a <- 10; b <- 3; c <- 4 f <- function(x) {exp(-a*x^3-b*x^2-c*x)} integrate(f,0,Inf) # works fine My difficulties start when I want to vectorize. # attempts to vectorize fail a <- seq(from=0,to=1,by=0.5) b <- seq(from=5,to=10,by=1) m <- seq(from=10,to=20,by=5) f2 <- function(x,a,b,c) {exp(-a*x^3-b*x^2-m*x)} fv <- Vectorize(integrate(f2,0,Inf),vectorize.args=c("a","b","m"),SIMPLIFY=TRUE) I want the result as a 3-d array with dimensions of the lengths of a, b and c. I have tried several variants but am not having much luck. Will appreciate any help that I can get. Thanks,Ravi Sutradhara         [[alternative HTML version deleted]] ______________________________________________ [hidden email] mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Open this post in threaded view
|

Re: vectorizing the integrate function

 Ravi: First of all, you're calling Vectorize incorrectly. The first argument must be a function *name*, not a function call. Here's what you need to do (and thanks for the reprex -- wouldn't have been able to help without it!): f2 <- function(a,b,c)integrate(function(x){exp(-a*x^3-b*x^2-c*x)},lower =0,upper = Inf) fv <- Vectorize(f2,vectorize.args=c("a","b","c"),SIMPLIFY=TRUE) Second of all, as you may realize, the vectorization uses mapply() and so vectorizes the c(a,b,c) triplet "in parallel," which will "recycle" shorter arguments to the length of longer. Hence you will get 6 results from your example: > fv(a,b,m)              [,1]         [,2]         [,3]        [,4]         [,5]   [,6] value        0.09207851   0.0635289    0.04837997  0.08856628   0.06224138   0.04777941 abs.error    3.365173e-08 3.108388e-06 1.00652e-09 3.284876e-09 1.796619e-08 6.348142e-09 subdivisions 2            1            2           2            2  2 message      "OK"         "OK"         "OK"        "OK"         "OK"   "OK" call         Expression   Expression   Expression  Expression   Expression   Expression Cheers, Bert Bert Gunter "The trouble with having an open mind is that people keep coming along and sticking things into it." -- Opus (aka Berkeley Breathed in his "Bloom County" comic strip ) On Sat, Aug 10, 2019 at 10:20 AM ravi via R-help <[hidden email]> wrote: > Hi all,I am having some difficulties in vectorizing the integrate > function. Let me explain with an example. > a <- 10; b <- 3; c <- 4 > f <- function(x) {exp(-a*x^3-b*x^2-c*x)} > integrate(f,0,Inf) # works fine > > My difficulties start when I want to vectorize. > > # attempts to vectorize fail > a <- seq(from=0,to=1,by=0.5) > b <- seq(from=5,to=10,by=1) > m <- seq(from=10,to=20,by=5) > f2 <- function(x,a,b,c) {exp(-a*x^3-b*x^2-m*x)} > fv <- > Vectorize(integrate(f2,0,Inf),vectorize.args=c("a","b","m"),SIMPLIFY=TRUE) > > I want the result as a 3-d array with dimensions of the lengths of a, b > and c. I have tried several variants but am not having much luck. Will > appreciate any help that I can get. > Thanks,Ravi Sutradhara > > > > >         [[alternative HTML version deleted]] > > ______________________________________________ > [hidden email] mailing list -- To UNSUBSCRIBE and more, see > 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. >         [[alternative HTML version deleted]] ______________________________________________ [hidden email] mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Open this post in threaded view
|

Re: vectorizing the integrate function

 Bert,Thanks a lot for your help. I have a few follow-up questions (shown in the comment lines).  Numerical values and error for a (i) vector and (ii) array. a <- seq(from=0,to=1,by=0.5) b <- seq(from=5,to=10,by=1) m <- seq(from=10,to=20,by=5) f2 <- function(a,b,m)integrate(function(x){exp(-a*x^3-b*x^2-m*x)},lower=0,upper = Inf) fv <- Vectorize(f2,vectorize.args=c("a","b","m"),SIMPLIFY=TRUE) r1 <- fv(a,b,m) # How to get just the numerical results as an array # for scalar values of a,b and m, I can get it with #s <- integrate(f,0,Inf)\$value # How do I get this for a vector # Finally, if I want an array for all the values of a, b and m, how should I proceed? r <- array(0,c(3,6,3)) r.error <- array(0,c(3,6,3)) # I want the results of the vectorized integrate function in the above arrays. How do I extract these values? Thanks,Ravi     On Saturday, 10 August 2019, 20:03:22 CEST, Bert Gunter <[hidden email]> wrote:      Ravi: First of all, you're calling Vectorize incorrectly. The first argument must be a function *name*, not a function call. Here's what you need to do (and thanks for the reprex -- wouldn't have been able to help without it!): f2 <- function(a,b,c)integrate(function(x){exp(-a*x^3-b*x^2-c*x)},lower =0,upper = Inf) fv <- Vectorize(f2,vectorize.args=c("a","b","c"),SIMPLIFY=TRUE) Second of all, as you may realize, the vectorization uses mapply() and so vectorizes the c(a,b,c) triplet "in parallel," which will "recycle" shorter arguments to the length of longer. Hence you will get 6 results from your example: > fv(a,b,m)              [,1]         [,2]         [,3]        [,4]         [,5]         [,6]         value        0.09207851   0.0635289    0.04837997  0.08856628   0.06224138   0.04777941   abs.error    3.365173e-08 3.108388e-06 1.00652e-09 3.284876e-09 1.796619e-08 6.348142e-09 subdivisions 2            1            2           2            2            2           message      "OK"         "OK"         "OK"        "OK"         "OK"         "OK"         call         Expression   Expression   Expression  Expression   Expression   Expression  Cheers,Bert Bert Gunter "The trouble with having an open mind is that people keep coming along and sticking things into it." -- Opus (aka Berkeley Breathed in his "Bloom County" comic strip ) On Sat, Aug 10, 2019 at 10:20 AM ravi via R-help <[hidden email]> wrote: Hi all,I am having some difficulties in vectorizing the integrate function. Let me explain with an example. a <- 10; b <- 3; c <- 4 f <- function(x) {exp(-a*x^3-b*x^2-c*x)} integrate(f,0,Inf) # works fine My difficulties start when I want to vectorize. # attempts to vectorize fail a <- seq(from=0,to=1,by=0.5) b <- seq(from=5,to=10,by=1) m <- seq(from=10,to=20,by=5) f2 <- function(x,a,b,c) {exp(-a*x^3-b*x^2-m*x)} fv <- Vectorize(integrate(f2,0,Inf),vectorize.args=c("a","b","m"),SIMPLIFY=TRUE) I want the result as a 3-d array with dimensions of the lengths of a, b and c. I have tried several variants but am not having much luck. Will appreciate any help that I can get. Thanks,Ravi Sutradhara         [[alternative HTML version deleted]] ______________________________________________ [hidden email] mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.           [[alternative HTML version deleted]] ______________________________________________ [hidden email] mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Open this post in threaded view
|

Re: vectorizing the integrate function

Open this post in threaded view
|

Re: vectorizing the integrate function

Open this post in threaded view
|

Re: vectorizing the integrate function

 In reply to this post by Bert Gunter-2  Rui,Thanks for your help in getting the result with for loops. That is useful, but I have a complicated integral and I am interested in improving the performance with the benefit of vectorization. In the solution from Bert, I would like to know how I can extract the numerical vector from the function fv. That is, extract the fv\$value as a vector. I can then perhaps try to improvise and extend it to arrays.Thanks,Ravi     On Saturday, 10 August 2019, 20:03:22 CEST, Bert Gunter <[hidden email]> wrote:      Ravi: First of all, you're calling Vectorize incorrectly. The first argument must be a function *name*, not a function call. Here's what you need to do (and thanks for the reprex -- wouldn't have been able to help without it!): f2 <- function(a,b,c)integrate(function(x){exp(-a*x^3-b*x^2-c*x)},lower =0,upper = Inf) fv <- Vectorize(f2,vectorize.args=c("a","b","c"),SIMPLIFY=TRUE) Second of all, as you may realize, the vectorization uses mapply() and so vectorizes the c(a,b,c) triplet "in parallel," which will "recycle" shorter arguments to the length of longer. Hence you will get 6 results from your example: > fv(a,b,m)              [,1]         [,2]         [,3]        [,4]         [,5]         [,6]         value        0.09207851   0.0635289    0.04837997  0.08856628   0.06224138   0.04777941   abs.error    3.365173e-08 3.108388e-06 1.00652e-09 3.284876e-09 1.796619e-08 6.348142e-09 subdivisions 2            1            2           2            2            2           message      "OK"         "OK"         "OK"        "OK"         "OK"         "OK"         call         Expression   Expression   Expression  Expression   Expression   Expression  Cheers,Bert Bert Gunter "The trouble with having an open mind is that people keep coming along and sticking things into it." -- Opus (aka Berkeley Breathed in his "Bloom County" comic strip ) On Sat, Aug 10, 2019 at 10:20 AM ravi via R-help <[hidden email]> wrote: Hi all,I am having some difficulties in vectorizing the integrate function. Let me explain with an example. a <- 10; b <- 3; c <- 4 f <- function(x) {exp(-a*x^3-b*x^2-c*x)} integrate(f,0,Inf) # works fine My difficulties start when I want to vectorize. # attempts to vectorize fail a <- seq(from=0,to=1,by=0.5) b <- seq(from=5,to=10,by=1) m <- seq(from=10,to=20,by=5) f2 <- function(x,a,b,c) {exp(-a*x^3-b*x^2-m*x)} fv <- Vectorize(integrate(f2,0,Inf),vectorize.args=c("a","b","m"),SIMPLIFY=TRUE) I want the result as a 3-d array with dimensions of the lengths of a, b and c. I have tried several variants but am not having much luck. Will appreciate any help that I can get. Thanks,Ravi Sutradhara         [[alternative HTML version deleted]] ______________________________________________ [hidden email] mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.           [[alternative HTML version deleted]] ______________________________________________ [hidden email] mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Open this post in threaded view
|

Re: vectorizing the integrate function

Open this post in threaded view
|

Re: vectorizing the integrate function

 In reply to this post by R help mailing list-2 Ravi: I believe you are under a misconception. The Vectorize() function vectorizes the function call, **but it does not vectorize the computation by moving loops down to the C level**, which is typically what is meant when it is recommended that R users use inbuilt vectorized functions when possible to move looping from the interpreted (R) level to the compiled (C level) code. That is, mapply() -- and thus Vectorize, which is just a wrapper for mapply -- is still looping at the interpreted level and therefore is essentially the same as explicit iterative or *apply loops. It may be a bit faster or slower, depending, but not the orders of magnitude faster that one can get with true vectorized (at the C level) code. Finally, if you want to just add the integration results to the allvals list, here's a slightly different one-liner to do it: > allvals\$vals <- unlist(with(allvals,fv(a,b,m))[1,]) > head(allvals)     a b  m       vals 1 0.0 5 10 0.09207851 2 0.5 5 10 0.09193111 3 1.0 5 10 0.09178672 4 0.0 6 10 0.09083412 5 0.5 6 10 0.09070066 6 1.0 6 10 0.09056960 or the one liner using the previous do.call() construction that just extracts the vector of values: > vals <- unlist(do.call(fv,allvals)[1,]) Cheers, Bert Bert Gunter "The trouble with having an open mind is that people keep coming along and sticking things into it." -- Opus (aka Berkeley Breathed in his "Bloom County" comic strip ) On Sat, Aug 10, 2019 at 3:56 PM ravi <[hidden email]> wrote: > Rui, > Thanks for your help in getting the result with for loops. That is useful, > but I have a complicated integral and I am interested in improving the > performance with the benefit of vectorization. > > In the solution from Bert, I would like to know how I can extract the > numerical vector from the function fv. That is, extract the fv\$value as a > vector. I can then perhaps try to improvise and extend it to arrays. > Thanks, > Ravi > > On Saturday, 10 August 2019, 20:03:22 CEST, Bert Gunter < > [hidden email]> wrote: > > > Ravi: > > First of all, you're calling Vectorize incorrectly. The first argument > must be a function *name*, not a function call. Here's what you need to do > (and thanks for the reprex -- wouldn't have been able to help without it!): > > f2 <- function(a,b,c)integrate(function(x){exp(-a*x^3-b*x^2-c*x)},lower > =0,upper = Inf) > fv <- Vectorize(f2,vectorize.args=c("a","b","c"),SIMPLIFY=TRUE) > > Second of all, as you may realize, the vectorization uses mapply() and so > vectorizes the c(a,b,c) triplet "in parallel," which will "recycle" shorter > arguments to the length of longer. Hence you will get 6 results from your > example: > > > fv(a,b,m) >              [,1]         [,2]         [,3]        [,4]         [,5] >   [,6] > value        0.09207851   0.0635289    0.04837997  0.08856628   0.06224138 >   0.04777941 > abs.error    3.365173e-08 3.108388e-06 1.00652e-09 3.284876e-09 > 1.796619e-08 6.348142e-09 > subdivisions 2            1            2           2            2 >    2 > message      "OK"         "OK"         "OK"        "OK"         "OK" >   "OK" > call         Expression   Expression   Expression  Expression   Expression >   Expression > > Cheers, > Bert > > Bert Gunter > > "The trouble with having an open mind is that people keep coming along and > sticking things into it." > -- Opus (aka Berkeley Breathed in his "Bloom County" comic strip ) > > > On Sat, Aug 10, 2019 at 10:20 AM ravi via R-help <[hidden email]> > wrote: > > Hi all,I am having some difficulties in vectorizing the integrate > function. Let me explain with an example. > a <- 10; b <- 3; c <- 4 > f <- function(x) {exp(-a*x^3-b*x^2-c*x)} > integrate(f,0,Inf) # works fine > > My difficulties start when I want to vectorize. > > # attempts to vectorize fail > a <- seq(from=0,to=1,by=0.5) > b <- seq(from=5,to=10,by=1) > m <- seq(from=10,to=20,by=5) > f2 <- function(x,a,b,c) {exp(-a*x^3-b*x^2-m*x)} > fv <- > Vectorize(integrate(f2,0,Inf),vectorize.args=c("a","b","m"),SIMPLIFY=TRUE) > > I want the result as a 3-d array with dimensions of the lengths of a, b > and c. I have tried several variants but am not having much luck. Will > appreciate any help that I can get. > Thanks,Ravi Sutradhara > > > > > >         [[alternative HTML version deleted]] > > ______________________________________________ > [hidden email] mailing list -- To UNSUBSCRIBE and more, see > 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. > >         [[alternative HTML version deleted]] ______________________________________________ [hidden email] mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.