As my post showed, it is a scaling issue. The function has so small a
peak that it is effectively 0 -- when scaled sensibly integrate then works
out of the box.
On Thu, 21 Aug 2008, Thomas Lumley wrote:
On Thu, 21 Aug 2008, Moshe Olshansky wrote:
The phenomenon is most likely caused by numerical errors. I do not know how
'integrate' works but numerical integration over a very long interval does
not look a good idea to me.
Numerical integration over a long (or infinite) interval is fine. The problem
here is that integrate appears not to find the peak of the integrand. The
integrand is very flat near zero, which makes things harder.
It seems to work well to split the integral somewhere not too far to the left
of its peak. Splits at 20, 50,100,or 200 give good agreement.
integrate(f,0,100)
0.0001600355 with absolute error < 6.5e-10
integrate(f,100,Inf)
5.355248e-06 with absolute error < 3.4e-06
integrate(f,0,200)
0.0001653797 with absolute error < 3.1e-05
integrate(f,200,Inf)
3.141137e-13 with absolute error < 1.2e-13
integrate(f,0,50)
2.702318e-05 with absolute error < 1.3e-16
integrate(f,50, Inf)
0.0001383181 with absolute error < 8e-05
Transforming the integral is only going to be helpful if you have a good idea
of what it looks like and what the integrate() function is expecting. An
automated transformation to a bounded interval won't help since that is the
first thing the underlying Fortran code does.
-thomas
I would do the following:
f1<-function(x){
return(dchisq(x,9,77)*((13.5/x)^5)*exp(-13.5/x))
}
f2<-function(y){
return(dchisq(1/y,9,77)*((13.5*y)^5)*exp(-13.5*y)/y^2)
}
# substitute y = 1/x
I1 <- integrate(f1,0,1,abs.tol=abs.tol=1.0e-20*.Machine$double.eps^0.25)
I2 <- integrate(f2,0,1,abs.tol=abs.tol=1.0e-20*.Machine$double.eps^0.25)
--- On Thu, 21/8/08, A <[EMAIL PROTECTED]> wrote:
From: A <[EMAIL PROTECTED]>
Subject: [R] Help Regarding 'integrate'
To: r-help@r-project.org
Received: Thursday, 21 August, 2008, 4:44 PM
I have an R function defined as:
f<-function(x){
return(dchisq(x,9,77)*((13.5/x)^5)*exp(-13.5/x))
}
Numerically integrating this function, I observed a couple
of things:
A) Integrating the function f over the entire positive real
line gives an
error:
integrate(f,0,Inf)
Error in integrate(f, 0, Inf) : the integral is probably
divergent
B) Increasing the interval of integration actually
decreases the value of
the integral:
integrate(f,0,10^5)
9.813968e-06 with absolute error < 1.9e-05
integrate(f,0,10^6)
4.62233e-319 with absolute error < 4.6e-319
Since the function f is uniformly positive, B) can not
occur. Also, the
theory tells me that the infinite integral actually exists
and is finite, so
A) can not occur. That means there are certain problems
with the usage of
function 'integrate' which I do not understand. The
help document tells me
that 'integrate' uses quadrature approximation to
evaluate integrals
numerically. Since I do not come from the numerical methods
community, I
would not know the pros and cons of various methods of
quadrature
approximation. One naive way that I thought of evaluating
the above integral
was by first locating the maximum of the function (may be
by using
'optimize' in R) and then applying the
Simpson's rule to an interval around
the maximum. However, I am sure that the people behind the
R project and
other users have much better ideas, and I am sure the
'correct' method has
already been implemented in R. Therefore, I would
appreciate if someone can
help me find it.
Thanks
[[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.
Thomas Lumley Assoc. Professor, Biostatistics
[EMAIL PROTECTED] University of Washington, Seattle
______________________________________________
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.
--
Brian D. Ripley, [EMAIL PROTECTED]
Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/
University of Oxford, Tel: +44 1865 272861 (self)
1 South Parks Road, +44 1865 272866 (PA)
Oxford OX1 3TG, UK Fax: +44 1865 272595
______________________________________________
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.
