Re: [R] Integration in R

[David Winsemius]

On Jan 8, 2013, at 1:31 PM, David Winsemius wrote:



On Jan 8, 2013, at 1:07 PM, Berend Hasselman wrote:



On 08-01-2013, at 22:00, Berend Hasselman  wrote:


…...
David implemented the condition by multiplying by x[1]in a numeric context is 0 when x[1]=x[2].


OOPS!! Reverse the 0 and the 1 in that sentence (TRUE becomes 1 and  
FALSE becomes 0)


I didn't worry about which direction was used for the inequality in  
this case because of consideration of symmetry. (Thinking about it  
now, I still think I managed to choose the correct direction.)   
Obviously such a cavalier attitude should not be carried over to the  
general situation.


I hit the send button by mistake. (I originally misread Berend's  
correction.) I'm only posting this to make clear that it was myself,  
and not Berend, that I was accusing of having been cavalier.


--

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.

Re: [R] Integration in R

[David Winsemius]

On Jan 8, 2013, at 1:07 PM, Berend Hasselman wrote:



On 08-01-2013, at 22:00, Berend Hasselman  wrote:


…...
David implemented the condition by multiplying by x[1]in a numeric context is 0 when x[1]=x[2].


OOPS!! Reverse the 0 and the 1 in that sentence (TRUE becomes 1 and  
FALSE becomes 0)


I didn't worry about which direction was used for the inequality in  
this case because of consideration of symmetry. (Thinking about it  
now, I still think I managed to choose the correct direction.)   
Obviously such a cavalier attitude should not be carried over to the  
general situation.


--

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.

Re: [R] Integration in R

[Berend Hasselman]

On 08-01-2013, at 22:00, Berend Hasselman  wrote:

> …...
> David implemented the condition by multiplying by x[1] numeric context is 0 when x[1]=x[2].

OOPS!! Reverse the 0 and the 1 in that sentence (TRUE becomes 1 and FALSE 
becomes 0)

Berend

__
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.

Re: [R] Integration in R

[Berend Hasselman]

On 08-01-2013, at 19:51, Naser Jamil  wrote:

> Thanks. But then how to implement condition like 0 happy to know that.
> 

David implemented the condition by multiplying by x[1]=x[2]. That is what your requirement 
does.
The condition 0x2 which is what your inequality implies.


Berend

> On 8 January 2013 18:41, David Winsemius  wrote:
> 
>> Please reply on list.
>> 
>> 
>> On Jan 8, 2013, at 10:27 AM, Naser Jamil wrote:
>> 
>> Hi David,
>>> x[2] is the second variable, x2. It comes from the condition 0>> 
>> 
>> No, it doesn't come from those conditions. It is being grabbed from some
>> "x"-named object that exists in your workspace.
>> 
>> If your limits were 7 in both dimensions, then the code should be:
>> 
>> adaptIntegrate(f, lowerLimit = c(0, 0), upperLimit = c(7,7))
>> #
>> $integral
>> [1] 228.6667
>> 
>> (At this point I was trusting R's calculus abilities more than yours. I
>> wasn't too trusting of mine either, and so tried seeing if Wolfram Alpha
>> would accept this expression:
>> integrate 2/3 (x+y) over 0< x<7,  0> 
>> ; which it did and calculating the decimal expansion of the exact fraction:
>> 
>>> 686/3
>> [1] 228.6667
>>> 
>> 
>> --
>> David.
>> 
>> 
>> 
>> 
>>> Thanks.
>>> 
>>> On 8 January 2013 18:11, David Winsemius  wrote:
>>> 
>>> On Jan 8, 2013, at 9:43 AM, Naser Jamil wrote:
>>> 
>>> Hi R-users.
>>> 
>>> I'm having difficulty with an integration in R via
>>> the package "cubature". I'm putting it with a simple example here.  I wish
>>> to integrate a function like:
>>> f(x1,x2)=2/3*(x1+x2) in the interval 0>> by hand and got 114.33, but the following R code is giving me 102.6667.
>>> 
>>> --**--**---
>>> library(cubature)
>>> f<-function(x) { 2/3 * (x[1] + x[2] ) }
>>> adaptIntegrate(f, lowerLimit = c(0, 0), upperLimit = c(x[2],7))
>>> 
>>> 
>>> What is x[2]?  On my machine it was 0.0761, so I obviously got a
>>> different answer.
>>> 
>>> --
>>> David Winsemius, MD
>>> Alameda, CA, USA
>>> 
>>> 
>>> 
>> David Winsemius, MD
>> Alameda, CA, USA
>> 
>> 
> 
>   [[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.

__
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.

Re: [R] Integration in R

[Naser Jamil]

Thanks. But then how to implement condition like 0 wrote:

> Please reply on list.
>
>
> On Jan 8, 2013, at 10:27 AM, Naser Jamil wrote:
>
>  Hi David,
>> x[2] is the second variable, x2. It comes from the condition 0>
>
> No, it doesn't come from those conditions. It is being grabbed from some
> "x"-named object that exists in your workspace.
>
> If your limits were 7 in both dimensions, then the code should be:
>
> adaptIntegrate(f, lowerLimit = c(0, 0), upperLimit = c(7,7))
> #
> $integral
> [1] 228.6667
>
> (At this point I was trusting R's calculus abilities more than yours. I
> wasn't too trusting of mine either, and so tried seeing if Wolfram Alpha
> would accept this expression:
>  integrate 2/3 (x+y) over 0< x<7,  0
> ; which it did and calculating the decimal expansion of the exact fraction:
>
> > 686/3
> [1] 228.6667
> >
>
> --
> David.
>
>
>
>
>> Thanks.
>>
>> On 8 January 2013 18:11, David Winsemius  wrote:
>>
>> On Jan 8, 2013, at 9:43 AM, Naser Jamil wrote:
>>
>> Hi R-users.
>>
>> I'm having difficulty with an integration in R via
>> the package "cubature". I'm putting it with a simple example here.  I wish
>> to integrate a function like:
>> f(x1,x2)=2/3*(x1+x2) in the interval 0> by hand and got 114.33, but the following R code is giving me 102.6667.
>>
>> --**--**---
>> library(cubature)
>> f<-function(x) { 2/3 * (x[1] + x[2] ) }
>> adaptIntegrate(f, lowerLimit = c(0, 0), upperLimit = c(x[2],7))
>>
>>
>> What is x[2]?  On my machine it was 0.0761, so I obviously got a
>> different answer.
>>
>> --
>> David Winsemius, MD
>> Alameda, CA, USA
>>
>>
>>
> David Winsemius, MD
> Alameda, CA, USA
>
>

[[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.

Re: [R] Integration in R

[David Winsemius]

On Jan 8, 2013, at 10:51 AM, Naser Jamil wrote:

Thanks. But then how to implement condition like 0be happy to know that.


Multiply the function by the conditional expression:

> f<-function(x) { 2/3 * (x[1] + x[2] )*(x[1] < x[2]) }
> adaptIntegrate(f, lowerLimit = c(0, 0), upperLimit = c(7,7))
$integral
[1] 114.

--
David.


On 8 January 2013 18:41, David Winsemius   
wrote:

Please reply on list.


On Jan 8, 2013, at 10:27 AM, Naser Jamil wrote:

Hi David,
x[2] is the second variable, x2. It comes from the condition  
0

No, it doesn't come from those conditions. It is being grabbed from  
some "x"-named object that exists in your workspace.


If your limits were 7 in both dimensions, then the code should be:

adaptIntegrate(f, lowerLimit = c(0, 0), upperLimit = c(7,7))
#
$integral
[1] 228.6667

(At this point I was trusting R's calculus abilities more than  
yours. I wasn't too trusting of mine either, and so tried seeing if  
Wolfram Alpha would accept this expression:

 integrate 2/3 (x+y) over 0< x<7,  0; which it did and calculating the decimal expansion of the exact  
fraction:


> 686/3
[1] 228.6667
>

--
David.




Thanks.

On 8 January 2013 18:11, David Winsemius   
wrote:


On Jan 8, 2013, at 9:43 AM, Naser Jamil wrote:

Hi R-users.

I'm having difficulty with an integration in R via
the package "cubature". I'm putting it with a simple example here.   
I wish

to integrate a function like:
f(x1,x2)=2/3*(x1+x2) in the interval 0by hand and got 114.33, but the following R code is giving me  
102.6667.


---
library(cubature)
f<-function(x) { 2/3 * (x[1] + x[2] ) }
adaptIntegrate(f, lowerLimit = c(0, 0), upperLimit = c(x[2],7))


What is x[2]?  On my machine it was 0.0761, so I obviously got a  
different answer.


--
David Winsemius, MD
Alameda, CA, USA



David Winsemius, MD
Alameda, CA, USA




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.

Re: [R] Integration in R

[David Winsemius]

Please reply on list.

On Jan 8, 2013, at 10:27 AM, Naser Jamil wrote:


Hi David,
x[2] is the second variable, x2. It comes from the condition  
0

No, it doesn't come from those conditions. It is being grabbed from  
some "x"-named object that exists in your workspace.


If your limits were 7 in both dimensions, then the code should be:

adaptIntegrate(f, lowerLimit = c(0, 0), upperLimit = c(7,7))
#
$integral
[1] 228.6667

(At this point I was trusting R's calculus abilities more than yours.  
I wasn't too trusting of mine either, and so tried seeing if Wolfram  
Alpha would accept this expression:

 integrate 2/3 (x+y) over 0< x<7,  0; which it did and calculating the decimal expansion of the exact  
fraction:


> 686/3
[1] 228.6667
>

--
David.




Thanks.

On 8 January 2013 18:11, David Winsemius   
wrote:


On Jan 8, 2013, at 9:43 AM, Naser Jamil wrote:

Hi R-users.

I'm having difficulty with an integration in R via
the package "cubature". I'm putting it with a simple example here.   
I wish

to integrate a function like:
f(x1,x2)=2/3*(x1+x2) in the interval 0by hand and got 114.33, but the following R code is giving me  
102.6667.


---
library(cubature)
f<-function(x) { 2/3 * (x[1] + x[2] ) }
adaptIntegrate(f, lowerLimit = c(0, 0), upperLimit = c(x[2],7))


What is x[2]?  On my machine it was 0.0761, so I obviously got a  
different answer.


--
David Winsemius, MD
Alameda, CA, USA




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.

Re: [R] Integration in R

[David Winsemius]

On Jan 8, 2013, at 9:43 AM, Naser Jamil wrote:


Hi R-users.

I'm having difficulty with an integration in R via
the package "cubature". I'm putting it with a simple example here.   
I wish

to integrate a function like:
f(x1,x2)=2/3*(x1+x2) in the interval 0by hand and got 114.33, but the following R code is giving me  
102.6667.


---
library(cubature)
f<-function(x) { 2/3 * (x[1] + x[2] ) }
adaptIntegrate(f, lowerLimit = c(0, 0), upperLimit = c(x[2],7))



What is x[2]?  On my machine it was 0.0761, so I obviously got a  
different answer.


--
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] Integration in R

[Naser Jamil]

Hi R-users.

I'm having difficulty with an integration in R via
the package "cubature". I'm putting it with a simple example here.  I wish
to integrate a function like:
f(x1,x2)=2/3*(x1+x2) in the interval 0https://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.

Re: [R] Integration in R

[Ravi Varadhan]

Send us a reproducible R code that shows what you actually tried.

Ravi Varadhan, Ph.D.
Assistant Professor
The Center on Aging and Health
Division of Geriatric Medicine & Gerontology
Johns Hopkins University
rvarad...@jhmi.edu
410-502-2619


[[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.

Re: [R] Integration in R

[William Dunlap]

But if you use a smaller sdlog value then integrate does get it wrong because 
it does not find the delta-like function hidden somewhere between 0 and 
infinity.
  > integrate(function(x){dlnorm(x,meanlog=-0.3,sdlog=0.00041)},0,Inf)
  0 with absolute error < 0
  > integrate(function(x){dlnorm(x,meanlog=-0.3,sdlog=0.00041)},0,1)
  0 with absolute error < 0
Tell it where the bulk of the mass is and it works
  > integrate(function(x)dlnorm(x,-0.3,0.00041), 0.738, 0.743, 
subdivisions=10^3)
  1 with absolute error < 4.4e-05

Bill Dunlap
Spotfire, TIBCO Software
wdunlap tibco.com


> -Original Message-
> From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On 
> Behalf
> Of Rolf Turner
> Sent: Wednesday, November 21, 2012 12:57 PM
> To: Rehena Sultana
> Cc: r-help@r-project.org
> Subject: Re: [R] Integration in R
> 
> On 21/11/12 22:26, Rehena Sultana wrote:
> > Dear R - Experts,
> >
> > I am trying to integrate lognormal distribution (mu = -0.3 and sigma2 = 
> > 0.00041.. )
> based on the some hypothetical data. But I am getting 0 as the result. I have 
> checked
> that my R-code is correct as code is giving me result for some other data.
> >
> > As I understand, when I am integrating some pdf with in the range of (0, 
> > Inf), I should
> not get 0. How to handle this kind of problem?
> >
> > Please help. Looking forward for your reply.
> 
> Reproducible example?
> 
> When I do
> 
> integrate(function(x){dlnorm(x,meanlog=-0.3,sdlog=sqrt(0.00041))},0,Inf)
> 
> I get
> 
>  1 with absolute error < 3.5e-06
> 
> No problema.
> 
> Why do you want to integrate it anyhow?  You know the answer is 1.
> 
> Or if you want the integral from 0 to x for some x < infinity, just use
> plnorm().  (???)
> 
>  cheers,
> 
>  Rolf Turner
> 
> __
> 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.

Re: [R] Integration in R

[Rolf Turner]

On 21/11/12 22:26, Rehena Sultana wrote:

Dear R - Experts,

I am trying to integrate lognormal distribution (mu = -0.3 and sigma2 = 
0.00041.. ) based on the some hypothetical data. But I am getting 0 as the 
result. I have checked that my R-code is correct as code is giving me result 
for some other data.

As I understand, when I am integrating some pdf with in the range of (0, Inf), 
I should not get 0. How to handle this kind of problem?

Please help. Looking forward for your reply.


Reproducible example?

When I do

integrate(function(x){dlnorm(x,meanlog=-0.3,sdlog=sqrt(0.00041))},0,Inf)

I get

1 with absolute error < 3.5e-06

No problema.

Why do you want to integrate it anyhow?  You know the answer is 1.

Or if you want the integral from 0 to x for some x < infinity, just use
plnorm().  (???)

cheers,

Rolf Turner

__
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] Integration in R

[Rehena Sultana]

Dear R - Experts, 

I am trying to integrate lognormal distribution (mu = -0.3 and sigma2 = 
0.00041.. ) based on the some hypothetical data. But I am getting 0 as the 
result. I have checked that my R-code is correct as code is giving me result 
for some other data.

As I understand, when I am integrating some pdf with in the range of (0, Inf), 
I should not get 0. How to handle this kind of problem? 

Please help. Looking forward for your reply. 

Regards,
rehena
[[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.

Re: [R] Integration in R

[Berend Hasselman]

On 02-10-2012, at 20:50, Dereje Bacha  wrote:

> Hi
> 
> I am facing a problem of restricting an intercept of systems of equations.  
> Y1=f(X1,X2,X3)
> Y2=(X1,X2,X4)
>  
> I want to restrict an intercept of equation 2 equal to coefficient of X2 of 
> equation 1. 
> 

Please do not hijack a thread/conversation.
Do not reply to a thread with a completely different subject.
Start a new thread for a new subject.

It is completely unclear what you want.

Berend

__
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.

Re: [R] Integration in R

[Dereje Bacha]

Hi

I am facing a problem of restricting an intercept of systems of equations.  
Y1=f(X1,X2,X3)
Y2=(X1,X2,X4)
 
I want to restrict an intercept of equation 2 equal to coefficient of X2 
of equation 1. 

Please help
Dereje  



 From: Naser Jamil 
To: r-help@r-project.org 
Sent: Tuesday, October 2, 2012 10:23 AM
Subject: [R] Integration in R

Dear R-users,
I am facing problem with integrating in R a likelihood function which is a
function of four parameters. It's giving me the result at the end but
taking more than half an hour to run. I'm wondering is there any other
efficient way deal with. The following is my code. I am ready to provide
any other description of my function if you need to move forward.



library(cubature)
dose<-c(2,3,5)
y0<-c(2,1,0)
y1<-c(1,1,1)
y2<-c(0,1,2)

lf<-function (x) {
    v<-1
    for (i in 1:length(dose)) {
        psi0<-1/((1+exp(x[1]+x[2]*dose[i]))*(1+exp(x[3]+x[4]*dose[i])))

psi1<-exp(x[1]+x[2]*dose[i])/((1+exp(x[1]+x[2]*dose[i]))*(1+exp(x[3]+x[4]*dose[i])))
       v<-v*(psi0^y0[i])*(psi1^y1[i])*((1-psi0-psi1)^y2[i])
        }
       return(v)
        }

adaptIntegrate(lf, lowerLimit = c(-20, 0,-20,0), upperLimit = c(0,10,0,10))

-
Thanks for your attention.


Kind regards,

Jamil.

    [[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.
[[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.

Re: [R] Integration in R

[William Dunlap]

Another nitpick: don't use return() in the last statement.
It isn't needed, it looks like some other language, and
dropping it saves 8% of the time for the uncompiled code
(the compiler seems to get of it).

Bill Dunlap
Spotfire, TIBCO Software
wdunlap tibco.com


> -Original Message-
> From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On 
> Behalf
> Of Berend Hasselman
> Sent: Tuesday, October 02, 2012 11:19 AM
> To: Rui Barradas
> Cc: r-help@r-project.org; Naser Jamil
> Subject: Re: [R] Integration in R
> 
> 
> On 02-10-2012, at 20:01, Rui Barradas  wrote:
> 
> > Hello,
> >
> > Yes, it's possible to remove the loop. Since the loop is used to compute a 
> > running
> product and all we want is the final result, use the vectorized behavior of R 
> and a final
> ?prod().
> > Seedup: another 2x. And 4x2 == 8 == 1 [decimal] order of magnitude.
> >
> >
> > lf2 <-function (x) {
> >   v<-1
> >   x1 <- x[1]
> >   x2 <- x[2]
> >   x3 <- x[3]
> >   x4 <- x[4]
> >   z1 <- exp(x1+x2*dose)
> >   z2 <- exp(x3+x4*dose)
> >   psi0<-1/((1+z1)*(1+z2))
> >   psi1<-z1*psi0
> >   v <- (psi0^y0)*(psi1^y1)*((1-psi0-psi1)^y2)
> >   return( prod(v) )
> > }
> >
> > lf2.c <- cmpfun(lf2)
> >
> > Hope this helps,
> 
> 
> Wonderful. It certainly does help.
> A single nitpick: the v <- 1 at the start of the function can now be removed.
> 
> I got a speedup of 7.5 compared to the very first version lf1.
> 
> Berend
> 
> __
> 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.

Re: [R] Integration in R

[Rui Barradas]

Hello,

Em 02-10-2012 19:18, Berend Hasselman escreveu:

On 02-10-2012, at 20:01, Rui Barradas  wrote:


Hello,

Yes, it's possible to remove the loop. Since the loop is used to compute a 
running product and all we want is the final result, use the vectorized 
behavior of R and a final ?prod().
Seedup: another 2x. And 4x2 == 8 == 1 [decimal] order of magnitude.


lf2 <-function (x) {
   v<-1
   x1 <- x[1]
   x2 <- x[2]
   x3 <- x[3]
   x4 <- x[4]
   z1 <- exp(x1+x2*dose)
   z2 <- exp(x3+x4*dose)
   psi0<-1/((1+z1)*(1+z2))
   psi1<-z1*psi0
   v <- (psi0^y0)*(psi1^y1)*((1-psi0-psi1)^y2)
   return( prod(v) )
}

lf2.c <- cmpfun(lf2)

Hope this helps,


Wonderful. It certainly does help.
A single nitpick: the v <- 1 at the start of the function can now be removed.

Yes, I thought about removing it but in the end forgot to.


I got a speedup of 7.5 compared to the very first version lf1.

My system is a Windows 7, R 2.15.1.

Rui Barradas

sessionInfo()
R version 2.15.1 (2012-06-22)
Platform: i386-pc-mingw32/i386 (32-bit)

locale:
[1] LC_COLLATE=Portuguese_Portugal.1252 LC_CTYPE=Portuguese_Portugal.1252
[3] LC_MONETARY=Portuguese_Portugal.1252 LC_NUMERIC=C
[5] LC_TIME=Portuguese_Portugal.1252

attached base packages:
[1] compiler  stats graphics  grDevices utils datasets methods
[8] base

other attached packages:
[1] rbenchmark_1.0.0 cubature_1.1-1

loaded via a namespace (and not attached):
[1] fortunes_1.5-0 tools_2.15.1


Berend



__
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.

Re: [R] Integration in R

[Berend Hasselman]

On 02-10-2012, at 20:01, Rui Barradas  wrote:

> Hello,
> 
> Yes, it's possible to remove the loop. Since the loop is used to compute a 
> running product and all we want is the final result, use the vectorized 
> behavior of R and a final ?prod().
> Seedup: another 2x. And 4x2 == 8 == 1 [decimal] order of magnitude.
> 
> 
> lf2 <-function (x) {
>   v<-1
>   x1 <- x[1]
>   x2 <- x[2]
>   x3 <- x[3]
>   x4 <- x[4]
>   z1 <- exp(x1+x2*dose)
>   z2 <- exp(x3+x4*dose)
>   psi0<-1/((1+z1)*(1+z2))
>   psi1<-z1*psi0
>   v <- (psi0^y0)*(psi1^y1)*((1-psi0-psi1)^y2)
>   return( prod(v) )
> }
> 
> lf2.c <- cmpfun(lf2)
> 
> Hope this helps,


Wonderful. It certainly does help.
A single nitpick: the v <- 1 at the start of the function can now be removed.

I got a speedup of 7.5 compared to the very first version lf1.

Berend

__
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.

Re: [R] Integration in R

[Rui Barradas]

Hello,

Yes, it's possible to remove the loop. Since the loop is used to compute 
a running product and all we want is the final result, use the 
vectorized behavior of R and a final ?prod().

Seedup: another 2x. And 4x2 == 8 == 1 [decimal] order of magnitude.


lf2 <-function (x) {
   v<-1
   x1 <- x[1]
   x2 <- x[2]
   x3 <- x[3]
   x4 <- x[4]
   z1 <- exp(x1+x2*dose)
   z2 <- exp(x3+x4*dose)
   psi0<-1/((1+z1)*(1+z2))
   psi1<-z1*psi0
   v <- (psi0^y0)*(psi1^y1)*((1-psi0-psi1)^y2)
   return( prod(v) )
}

lf2.c <- cmpfun(lf2)

Hope this helps,

Rui Barradas
Em 02-10-2012 18:21, Berend Hasselman escreveu:

On 02-10-2012, at 17:23, Naser Jamil  wrote:


Dear R-users,
I am facing problem with integrating in R a likelihood function which is a
function of four parameters. It's giving me the result at the end but
taking more than half an hour to run. I'm wondering is there any other
efficient way deal with. The following is my code. I am ready to provide
any other description of my function if you need to move forward.



library(cubature)
dose<-c(2,3,5)
y0<-c(2,1,0)
y1<-c(1,1,1)
y2<-c(0,1,2)

lf<-function (x) {
v<-1
for (i in 1:length(dose)) {
psi0<-1/((1+exp(x[1]+x[2]*dose[i]))*(1+exp(x[3]+x[4]*dose[i])))

psi1<-exp(x[1]+x[2]*dose[i])/((1+exp(x[1]+x[2]*dose[i]))*(1+exp(x[3]+x[4]*dose[i])))
   v<-v*(psi0^y0[i])*(psi1^y1[i])*((1-psi0-psi1)^y2[i])
}
   return(v)
}

adaptIntegrate(lf, lowerLimit = c(-20, 0,-20,0), upperLimit = c(0,10,0,10))


There are several things you could do.

1. Use the compiler package to make a compiled version of your function.
2. rewrite your function by putting x[1] etc. in separate variables x1, x2,... 
so avoiding the [..] indexing. You can do the same for dose[i].
And also compiling this version of the function.
3. do not recompute expressions such as exp(x1+x2*dose.i) several times. Store 
the result in a temporary variable and use that variable.

With these functions

library(compiler)

lf.c <- cmpfun(lf)

lf1 <-function (x) {
v<-1
x1 <- x[1]
x2 <- x[2]
x3 <- x[3]
x4 <- x[4]
for (i in 1:length(dose)) {
dose.i <- dose[i]
z1 <- exp(x1+x2*dose.i)
z2 <- exp(x3+x4*dose.i)
psi0<-1/((1+z1)*(1+z2))
psi1<-z1*psi0
v<-v*(psi0^y0[i])*(psi1^y1[i])*((1-psi0-psi1)^y2[i])
}
return(v)
}

lf1.c <- cmpfun(lf1)

I tested relative speeds with this code (small tolerance and max. function 
evaluations)

library(rbenchmark)

f1 <- function() adaptIntegrate(lf   , lowerLimit = c(-20, 0,-20,0), upperLimit 
= c(0,10,0,10), tol=1e-3,maxEval=5)
f2 <- function() adaptIntegrate(lf.c , lowerLimit = c(-20, 0,-20,0), upperLimit 
= c(0,10,0,10), tol=1e-3,maxEval=5)
f3 <- function() adaptIntegrate(lf1  , lowerLimit = c(-20, 0,-20,0), upperLimit 
= c(0,10,0,10), tol=1e-3,maxEval=5)
f4 <- function() adaptIntegrate(lf1.c, lowerLimit = c(-20, 0,-20,0), upperLimit 
= c(0,10,0,10), tol=1e-3,maxEval=5)
benchmark(z1 <- f1(),z2 <- f2(), z3 <- f3(),z4 <- f4(),replications=1)

Result:


benchmark(z1 <- f1(),z2 <- f2(), z3 <- f3(),z4 <- f4(),replications=1)

 test replications elapsed relative user.self sys.self user.child
1 z1 <- f1()1   3.1974.535 3.1770.008  0
2 z2 <- f2()1   1.2401.759 1.2350.003  0
3 z3 <- f3()1   2.1713.079 2.1670.002  0
4 z4 <- f4()1   0.7051.000 0.7000.004  0

As you can see you should be able to get at least a fourfold speedup by using 
the compiled version of the optimized function.
I would certainly  set tol and maxEval to something reasonable initially.
Checking that z1, z2, z3, and z4 are equal is left to you.

Finally, it may even be possible to eliminate the for loop in your function but 
I'll leave that for someone else.

Berend

__
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.

Re: [R] Integration in R

[Berend Hasselman]

On 02-10-2012, at 17:23, Naser Jamil  wrote:

> Dear R-users,
> I am facing problem with integrating in R a likelihood function which is a
> function of four parameters. It's giving me the result at the end but
> taking more than half an hour to run. I'm wondering is there any other
> efficient way deal with. The following is my code. I am ready to provide
> any other description of my function if you need to move forward.
> 
> 
> 
> library(cubature)
> dose<-c(2,3,5)
> y0<-c(2,1,0)
> y1<-c(1,1,1)
> y2<-c(0,1,2)
> 
> lf<-function (x) {
>v<-1
>for (i in 1:length(dose)) {
>psi0<-1/((1+exp(x[1]+x[2]*dose[i]))*(1+exp(x[3]+x[4]*dose[i])))
> 
> psi1<-exp(x[1]+x[2]*dose[i])/((1+exp(x[1]+x[2]*dose[i]))*(1+exp(x[3]+x[4]*dose[i])))
>   v<-v*(psi0^y0[i])*(psi1^y1[i])*((1-psi0-psi1)^y2[i])
>}
>   return(v)
>}
> 
> adaptIntegrate(lf, lowerLimit = c(-20, 0,-20,0), upperLimit = c(0,10,0,10))
> 

There are several things you could do.

1. Use the compiler package to make a compiled version of your function.
2. rewrite your function by putting x[1] etc. in separate variables x1, x2,... 
so avoiding the [..] indexing. You can do the same for dose[i].
And also compiling this version of the function.
3. do not recompute expressions such as exp(x1+x2*dose.i) several times. Store 
the result in a temporary variable and use that variable.

With these functions

library(compiler)

lf.c <- cmpfun(lf)

lf1 <-function (x) {
   v<-1 
   x1 <- x[1]
   x2 <- x[2]
   x3 <- x[3]
   x4 <- x[4]
   for (i in 1:length(dose)) { 
   dose.i <- dose[i]
   z1 <- exp(x1+x2*dose.i)
   z2 <- exp(x3+x4*dose.i)
   psi0<-1/((1+z1)*(1+z2))
   psi1<-z1*psi0
   v<-v*(psi0^y0[i])*(psi1^y1[i])*((1-psi0-psi1)^y2[i])
   }
   return(v)
}

lf1.c <- cmpfun(lf1)

I tested relative speeds with this code (small tolerance and max. function 
evaluations)

library(rbenchmark)

f1 <- function() adaptIntegrate(lf   , lowerLimit = c(-20, 0,-20,0), upperLimit 
= c(0,10,0,10), tol=1e-3,maxEval=5)
f2 <- function() adaptIntegrate(lf.c , lowerLimit = c(-20, 0,-20,0), upperLimit 
= c(0,10,0,10), tol=1e-3,maxEval=5)
f3 <- function() adaptIntegrate(lf1  , lowerLimit = c(-20, 0,-20,0), upperLimit 
= c(0,10,0,10), tol=1e-3,maxEval=5)
f4 <- function() adaptIntegrate(lf1.c, lowerLimit = c(-20, 0,-20,0), upperLimit 
= c(0,10,0,10), tol=1e-3,maxEval=5)
benchmark(z1 <- f1(),z2 <- f2(), z3 <- f3(),z4 <- f4(),replications=1)

Result:

> benchmark(z1 <- f1(),z2 <- f2(), z3 <- f3(),z4 <- f4(),replications=1)
test replications elapsed relative user.self sys.self user.child
1 z1 <- f1()1   3.1974.535 3.1770.008  0
2 z2 <- f2()1   1.2401.759 1.2350.003  0
3 z3 <- f3()1   2.1713.079 2.1670.002  0
4 z4 <- f4()1   0.7051.000 0.7000.004  0

As you can see you should be able to get at least a fourfold speedup by using 
the compiled version of the optimized function.
I would certainly  set tol and maxEval to something reasonable initially.
Checking that z1, z2, z3, and z4 are equal is left to you.

Finally, it may even be possible to eliminate the for loop in your function but 
I'll leave that for someone else.

Berend

__
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] Integration in R

[Naser Jamil]

Dear R-users,
I am facing problem with integrating in R a likelihood function which is a
function of four parameters. It's giving me the result at the end but
taking more than half an hour to run. I'm wondering is there any other
efficient way deal with. The following is my code. I am ready to provide
any other description of my function if you need to move forward.



library(cubature)
dose<-c(2,3,5)
y0<-c(2,1,0)
y1<-c(1,1,1)
y2<-c(0,1,2)

lf<-function (x) {
v<-1
for (i in 1:length(dose)) {
psi0<-1/((1+exp(x[1]+x[2]*dose[i]))*(1+exp(x[3]+x[4]*dose[i])))

psi1<-exp(x[1]+x[2]*dose[i])/((1+exp(x[1]+x[2]*dose[i]))*(1+exp(x[3]+x[4]*dose[i])))
   v<-v*(psi0^y0[i])*(psi1^y1[i])*((1-psi0-psi1)^y2[i])
}
   return(v)
}

adaptIntegrate(lf, lowerLimit = c(-20, 0,-20,0), upperLimit = c(0,10,0,10))

-
Thanks for your attention.


Kind regards,

Jamil.

[[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.

Re: [R] Integration in R

[Ravi Varadhan]

You are missing basic algebra skills!  

You had:

myfunc<- function(x) {0.25*(9*x^4 + 6*x^2 + 1)}

This should be:

myfunc<- function(x) {0.25*(9*x^4 - 6*x^2 + 1)}


Ravi.
---
Ravi Varadhan, Ph.D.
Assistant Professor,
Division of Geriatric Medicine and Gerontology School of Medicine Johns
Hopkins University

Ph. (410) 502-2619
email: rvarad...@jhmi.edu

-Original Message-
From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On
Behalf Of Alaios
Sent: Monday, January 10, 2011 9:13 AM
To: r-help@r-project.org
Subject: [R] Integration in R

Dear all,

It has been   ages since I studied integration in college. Right now I 
try to recover all this kind of knowledge and then try to understand how
 integration works. 



Thus I am doing some first 'experiments' and I would like to request your
help and comments.



I have the function:



p2<-function(x){0.5*(3*x^2-1)}

# I found the square of p2 by using some pencil and paper.

# The result is inside myfunc below



myfunc<- function(x) {0.25*(9*x^4+6*x^2+1)}



# Below I made R to find the square of p2

p2sq<-function(x) {p2(x) * p2(x)}



# Now I am trying to integrate both two function at the same interval. #Both
functions should denote the square of p2



integrate(p2sq,-1,+1)

# returns 0.4 with absolute error < 4.4e-15



integrate(myfunc,-1,+1)

# returns 2.4 with absolute error < 2.7e-14



if there is no error in my calculations could you please explain me why 
the two integrations return different results. Might be that I am 
missing something from the theory. I would like to thank you for your 
help



Regards

Alex

__
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] Integration in R

[Alaios]

Dear all,

It has been   ages since I studied integration in college. Right now I 
try to recover all this kind of knowledge and then try to understand how
 integration works. 



Thus I am doing some first 'experiments' and I would like to request your help 
and comments.



I have the function:



p2<-function(x){0.5*(3*x^2-1)}

# I found the square of p2 by using some pencil and paper.

# The result is inside myfunc below



myfunc<- function(x) {0.25*(9*x^4+6*x^2+1)}



# Below I made R to find the square of p2

p2sq<-function(x) {p2(x) * p2(x)}



# Now I am trying to integrate both two function at the same interval. #Both 
functions should denote the square of p2



integrate(p2sq,-1,+1)

# returns 0.4 with absolute error < 4.4e-15



integrate(myfunc,-1,+1)

# returns 2.4 with absolute error < 2.7e-14



if there is no error in my calculations could you please explain me why 
the two integrations return different results. Might be that I am 
missing something from the theory. I would like to thank you for your 
help



Regards

Alex

__
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.