[R] Integration in R
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 0x1x27. To be sure I tried it
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))
---
I guess something is wrong with the way I have assigned limits
in my codes. May I seek some advice from you.
Many thanks.
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
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 0x1x27. To be sure I tried it
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
__
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
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
0x1x27.
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 x7, 0y7
; 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 dwinsem...@comcast.net
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 0x1x27. To be sure I tried it
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
__
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
On Jan 8, 2013, at 10:51 AM, Naser Jamil wrote:
Thanks. But then how to implement condition like 0x1x27? I would
be 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 dwinsem...@comcast.net
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
0x1x27.
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 x7, 0y7
; 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 dwinsem...@comcast.net
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 0x1x27. To be sure I tried it
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
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
Thanks. But then how to implement condition like 0x1x27? I would be
happy to know that.
On 8 January 2013 18:41, David Winsemius dwinsem...@comcast.net 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 0x1x27.
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 x7, 0y7
; 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 dwinsem...@comcast.net 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 0x1x27. To be sure I tried it
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
On 08-01-2013, at 19:51, Naser Jamil jamilnase...@gmail.com wrote:
Thanks. But then how to implement condition like 0x1x27? I would be
happy to know that.
David implemented the condition by multiplying by x[1]x[2] which in a numeric
context is 0 when x[1]x[2] en 1 when x[1]=x[2]. That is what your requirement
does.
The condition 0x1x27 has been split into three separate parts
0x17
0x27
x1x2
The first two can be converted to arguments of adaptIntegrate.
The last is inserted into the function definition effectively making the
function return 0 when x1x2 which is what your inequality implies.
Berend
On 8 January 2013 18:41, David Winsemius dwinsem...@comcast.net 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 0x1x27.
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 x7, 0y7
; 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 dwinsem...@comcast.net 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 0x1x27. To be sure I tried it
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
On 08-01-2013, at 22:00, Berend Hasselman b...@xs4all.nl wrote: …... David implemented the condition by multiplying by x[1]x[2] which in a numeric context is 0 when x[1]x[2] en 1 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
On Jan 8, 2013, at 1:07 PM, Berend Hasselman wrote: On 08-01-2013, at 22:00, Berend Hasselman b...@xs4all.nl wrote: …... David implemented the condition by multiplying by x[1]x[2] which in a numeric context is 0 when x[1]x[2] en 1 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
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 b...@xs4all.nl wrote: …... David implemented the condition by multiplying by x[1]x[2] which in a numeric context is 0 when x[1]x[2] en 1 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.
[R] Integration in R
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
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.
Re: [R] Integration in R
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
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.edumailto: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.
[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.
Re: [R] Integration in R
On 02-10-2012, at 17:23, Naser Jamil jamilnase...@gmail.com 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.
Re: [R] Integration in R
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 jamilnase...@gmail.com 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
On 02-10-2012, at 20:01, Rui Barradas ruipbarra...@sapo.pt 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
Hello,
Em 02-10-2012 19:18, Berend Hasselman escreveu:
On 02-10-2012, at 20:01, Rui Barradas ruipbarra...@sapo.pt 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
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 ruipbarra...@sapo.pt 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
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 jamilnase...@gmail.com
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
On 02-10-2012, at 20:50, Dereje Bacha d_ba...@yahoo.com 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.
[R] Integration between R and Microsoft SQL Server?
Oracle has its own R integration called Oracle R Enterprise: http://www.oracle.com/technetwork/database/options/advanced-analytics/r-enterprise/index.html Is there a same kind of integration available for Microsoft SQL Server? If there is not, does MS have any plans to integrate R and SQL Server database engine? -J __ 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 between R and Microsoft SQL Server?
On 12-09-11 3:30 AM, johannes rara wrote: Oracle has its own R integration called Oracle R Enterprise: http://www.oracle.com/technetwork/database/options/advanced-analytics/r-enterprise/index.html Is there a same kind of integration available for Microsoft SQL Server? If there is not, does MS have any plans to integrate R and SQL Server database engine? I think you need to ask the second question on the Microsoft-SQL-Server-help mailing list (or whatever support forum MS uses), not R-help. Duncan Murdoch __ 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
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.
Re: [R] Integration in R
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.
