Re: [R] numerical integration
My goodness! Did you try ?integrate ?
Bert Gunter
Genentech Nonclinical Biostatistics
-Original Message-
From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On
Behalf Of Roslina Zakaria
Sent: Thursday, September 10, 2009 3:36 PM
To: r-help@r-project.org
Subject: [R] numerical integration
Hi r-users,
Can I do a numerical integration in R to solve for F(z)- integral_0^z {f(t)
dt} = 0 where F(z) is the CDF and f(t) is the pdf? What package can I use?
Thank you so much for any help given.
[[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] Numerical Integration Problems
Hi, You may want to try the double exponential transformation on the numerator and the denominator on this one. The method is described in detail here: http://projecteuclid.org/DPubS?service=UI&version=1.0&verb=Display&handle=euclid.prims/1145474600 If you want to give it a shot outside R there are a couple of C/C++ implementations out there: http://www.codeproject.com/KB/recipes/FastNumericalIntegration.aspx http://www.crbond.com/download/deint.c I have used the second one in a project of mine that is similar to the one you describe (calculation of posterior expectations in Bayesian inference), and I can say that it works reasonably well. As a side note, the DE method is one (actually the default) of the numerical integration methods implemented by Mathematica's NIntegrate function. It would be nice if someone found the time to augment R's integration facilities with this (and possibly other) numerical integration methods. Christos Argyropoulos University of Pittsburgh Medical Center _ s. It's easy! aspx&mkt=en-us __ 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] Numerical integration problem
Hi Marcus,
I always use a smaller error tolerance in `integrate' than the default
value. I generally use 1.e-07, whereas the default is only about 1.e-04.
Sometimes you may also need to increase the number of subdivisions from its
default value of 100.
Your problem disappears if you use a smaller tolerance for convergence:
> integrate(hazard, 0, Inf, v=0.1, gam=pi*231/200,
pos=list(r=5,th=atan(3/4)), rel.tol=1.e-07, a=10, b0=10, bt=0.1)
0.0002290673 with absolute error < 1.9e-11
>
Ravi.
---
Ravi Varadhan, Ph.D.
Assistant Professor, The Center on Aging and Health
Division of Geriatric Medicine and Gerontology
Johns Hopkins University
Ph: (410) 502-2619
Fax: (410) 614-9625
Email: rvarad...@jhmi.edu
Webpage:
http://www.jhsph.edu/agingandhealth/People/Faculty_personal_pages/Varadhan.h
tml
-Original Message-
From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On
Behalf Of Marcus Rowcliffe
Sent: Wednesday, September 23, 2009 11:31 AM
To: r-help@r-project.org
Subject: [R] Numerical integration problem
Hi there
I'm trying to construct a model of mortality risk in 2D space that
requires numerical integration of a hazard function, for which I'm using
the integrate function. I'm occasionally encountering parameter
combinations that cause integrate to terminate with error "Error in
integrate... the integral is probably divergent", which I'm not sure how
to interpret. The problem only crops up for a tiny part of the input
parameter space, but if I can't get this to work for the whole space I'm
a bit stuck! Plotting the integrand shows no sign of obvious problems
such as flatness, extreme spiking or non-finite area, and a tiny tweak
to the input is enough to get the integration to work, even though the
tweaked function looks essentially the same as the one that fails. Below
is some code that demonstrates the problem on my machine running R
2.9.0. Any suggestions on what might be causing this and what I can do
to avoid it would be very gratefully received.
Hoping for some insights
Marcus
##
#Problem example
#Works fine:
integrate(hazard, 0,Inf, v=0.1, gam=pi*232/200,
pos=list(r=5,th=atan(3/4)), a=10, b0=10, bt=0.1)
#Gives error "...integral probably divergent":
integrate(hazard, 0,Inf, v=0.1, gam=pi*231/200,
pos=list(r=5,th=atan(3/4)), a=10, b0=10, bt=0.1)
#Plot the integrands - doesn't look obviously problematic
h <-
hazard(seq(0,500,0.1),0.1,pi*231/200,pos=list(r=5,th=atan(3/4)),10,10,0.
1)
plot(seq(0,500,0.1),h,type="l",col=2)
h <-
hazard(seq(0,500,0.1),0.1,pi*232/200,pos=list(r=5,th=atan(3/4)),10,10,0.
1)
lines(seq(0,500,0.1),h)
#Functions used:
#hazard at a given point pos (a list of polar coordinates: distance r
and angle th); a, b0 and bt are model parameters
point.hazard <- function(pos,a,b0,bt) a * exp(-(pos$r^2/(2*b0^2))) *
exp(-(pos$th^2/(2*bt^2)))
#point.hazard for a point related to input point pos by time t and speed
v
hazard <- function(t, v, gam, pos, a, b0, bt)
{ pos2 <- zredef(pos,-t*v,gam)
pos2$th[pos2$th>pi] <- 2*pi-pos2$th[pos2$th>pi]
point.hazard(pos2,a,b0,bt)
}
#Returns a list of polar co-ordinates for a point defined by distance m
in direction gam from starting point pos
zredef <- function(pos, m, gam)
{ x <- pos$r*sin(pos$th) + m*sin(gam)
y <- pos$r*cos(pos$th) + m*cos(gam)
r <- (x^2 + y^2)^0.5
th <- atan(x/y)
th[x==0] <- 0
th[y<0] <- th[y<0] + pi
th[x<0 & y>0] <- th[x<0 & y>0] + 2*pi
list(r=r,th=th)
}
The Zoological Society of London is incorporated by Royal Charter
Principal Office England. Company Number RC000749
Registered address:
Regent's Park, London, England NW1 4RY
Registered Charity in England and Wales no. 208728
_
This e-mail has been sent in confidence to the named add...{{dropped:8}}
__
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] Numerical Integration in 1D
On Fri, 7 Mar 2008, Max wrote:
> Dear UseRs,
>
> I'm curious about the derivative of n!.
>
> We know that Gamma(n+1)=n! So when on takes the derivative of
> Gamma(n+1) we get Int(ln(x)*exp(-x)*x^n,x=0..Inf).
>
> I've tried code like
>
>> integrand<-function(x) {log(x)*exp(x)*x^n}
>> integrate(integrand,lower=0,upper=Inf)
>
> It seems that R doesn't like to integrate for any n, and I was
> wondering if anyone knew a way around this?
ln(x) e^x x^n is not integrable on (0, Inf). You presumably slipped over
a minus sign, but your definition of gamma(n) is wrong -- see ?gamma.
integrate(function(x) exp(-x)*x^n, lower=0, upper=Inf)
will work for gamma(n+1).
--
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, UKFax: +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.
Re: [R] Numerical Integration in 1D
Hi Max,
The analytic integral \int _0 ^\Inf exp(-t) t^n log(t) might not converge
because the integrand tends to -Inf as t -> 0.
So, here is a numerical approach to estimating the derivative of the gamma
function:
library(numDeriv)
fx <- function(x, n) exp(-x) * x^n
gf <- function(n) {integrate(fx, lower=0, upper=Inf, n=n)$val}
> grad(x=3, func=gf)
[1] 7.536706
>
> grad(x=10, func=gf)
[1] 8534040
>
Best,
Ravi.
---
Ravi Varadhan, Ph.D.
Assistant Professor, The Center on Aging and Health
Division of Geriatric Medicine and Gerontology
Johns Hopkins University
Ph: (410) 502-2619
Fax: (410) 614-9625
Email: [EMAIL PROTECTED]
Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html
-Original Message-
From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On
Behalf Of Max
Sent: Friday, March 07, 2008 1:41 PM
To: [EMAIL PROTECTED]
Subject: [R] Numerical Integration in 1D
Dear UseRs,
I'm curious about the derivative of n!.
We know that Gamma(n+1)=n! So when on takes the derivative of
Gamma(n+1) we get Int(ln(x)*exp(-x)*x^n,x=0..Inf).
I've tried code like
> integrand<-function(x) {log(x)*exp(x)*x^n}
> integrate(integrand,lower=0,upper=Inf)
It seems that R doesn't like to integrate for any n, and I was
wondering if anyone knew a way around this?
-Max
__
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] Numerical Integration in 1D
Prof Brian Ripley formulated on Friday :
> On Fri, 7 Mar 2008, Max wrote:
>
>> Dear UseRs,
>>
>> I'm curious about the derivative of n!.
>>
>> We know that Gamma(n+1)=n! So when on takes the derivative of
>> Gamma(n+1) we get Int(ln(x)*exp(-x)*x^n,x=0..Inf).
>>
>> I've tried code like
>>
>>> integrand<-function(x) {log(x)*exp(x)*x^n}
>>> integrate(integrand,lower=0,upper=Inf)
>>
>> It seems that R doesn't like to integrate for any n, and I was
>> wondering if anyone knew a way around this?
>
> ln(x) e^x x^n is not integrable on (0, Inf). You presumably slipped over
> a minus sign, but your definition of gamma(n) is wrong -- see ?gamma.
>
> integrate(function(x) exp(-x)*x^n, lower=0, upper=Inf)
>
> will work for gamma(n+1).
I did miss a minus sign in the integration, which explains part of my
problems. The function of interest is the derivative of Gamma(n+1) with
respect to n, but obviously integrated over x from 0 to Infinity.
What happens now is:
> integrand<-function(x) {log(x)*exp(-x)*x^n}
> integrate(integrand,lower=0,upper=Inf)
Error in f(x, ...) : object "n" not found
Any ideas on how to get around this error?
__
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] Numerical Integration in 1D
Hi max,
Prof. Ripley is right. Your problem is that you missed a (-) sign in the
exponential. Here is a demonstration showing the agreement between
numerical and analytical results:
gx <- function(x, n) exp(-x) * x^n * log(x)
df <- function(n) {integrate(gx, lower=0, upper=Inf, n=n)$val}
library(numDeriv)
fx <- function(x, n) exp(-x) * x^n
gf <- function(n) {integrate(fx, lower=0, upper=Inf, n=n)$val}
> grad(x=6, func=gf)
[1] 1348.405
>
> df(6)
[1] 1348.405
>
Ravi.
---
Ravi Varadhan, Ph.D.
Assistant Professor, The Center on Aging and Health
Division of Geriatric Medicine and Gerontology
Johns Hopkins University
Ph: (410) 502-2619
Fax: (410) 614-9625
Email: [EMAIL PROTECTED]
Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html
-Original Message-
From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On
Behalf Of Max
Sent: Friday, March 07, 2008 1:41 PM
To: [EMAIL PROTECTED]
Subject: [R] Numerical Integration in 1D
Dear UseRs,
I'm curious about the derivative of n!.
We know that Gamma(n+1)=n! So when on takes the derivative of
Gamma(n+1) we get Int(ln(x)*exp(-x)*x^n,x=0..Inf).
I've tried code like
> integrand<-function(x) {log(x)*exp(x)*x^n}
> integrate(integrand,lower=0,upper=Inf)
It seems that R doesn't like to integrate for any n, and I was
wondering if anyone knew a way around this?
-Max
__
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] Numerical Integration in 1D
On Fri, 7 Mar 2008, Max wrote:
> Prof Brian Ripley formulated on Friday :
>> On Fri, 7 Mar 2008, Max wrote:
>>
>>> Dear UseRs,
>>>
>>> I'm curious about the derivative of n!.
>>>
>>> We know that Gamma(n+1)=n! So when on takes the derivative of
>>> Gamma(n+1) we get Int(ln(x)*exp(-x)*x^n,x=0..Inf).
>>>
>>> I've tried code like
>>>
integrand<-function(x) {log(x)*exp(x)*x^n}
integrate(integrand,lower=0,upper=Inf)
>>>
>>> It seems that R doesn't like to integrate for any n, and I was
>>> wondering if anyone knew a way around this?
>>
>> ln(x) e^x x^n is not integrable on (0, Inf). You presumably slipped over
>> a minus sign, but your definition of gamma(n) is wrong -- see ?gamma.
>>
>> integrate(function(x) exp(-x)*x^n, lower=0, upper=Inf)
>>
>> will work for gamma(n+1).
>
> I did miss a minus sign in the integration, which explains part of my
> problems. The function of interest is the derivative of Gamma(n+1) with
> respect to n, but obviously integrated over x from 0 to Infinity.
And you said n!, so n must be integer and you cannot differentiate a
function of a integer argument.
If you are interested in the derivative of gamma(x), check out ?digamma.
> What happens now is:
>
>> integrand<-function(x) {log(x)*exp(-x)*x^n}
>> integrate(integrand,lower=0,upper=Inf)
> Error in f(x, ...) : object "n" not found
>
> Any ideas on how to get around this error?
Set 'n' in the evaluation environment.
E.g.
> n <- 3
> integrate(integrand, lower=0, upper=Inf)
7.536706 with absolute error < 4.7e-06
> digamma(4)*gamma(4)
[1] 7.536706
--
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, UKFax: +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.
Re: [R] numerical integration of a ftn of 2 variables
On Tue, Feb 19, 2008 at 11:07 PM, Chris Rhoads
<[EMAIL PROTECTED]> wrote:
> To start, let me confess to not being an experienced programmer, although I
> have used R fairly
> extensively in my work as a
> graduate student in statistics.
>
> I wish to find the root of a function of two variables that is defined by an
> integral which must be
> evaluated numerically.
>
> So the problem I want to solve is of the form: Find k such that f(k)=0,
> where f(y) = int_a^b
> g(x,y) dx. Again, the integral
> involved must be done numerically.
>
> I'm told by a friend who knows programming, but not R, that what I need to
> do is create something
> like a "local environment"
> within which I could create a placeholder for x. So I want to make
> something like the following work.
>
> f(var) <- function(var) {
>
> cons <- var
>
> g <- function(x,cons) {h(x,cons)}
>
> ret <- function(cons) integrate(g(x,cons),a,b)$value
> ret
> }
>
> I could then use (e.g.) a Newton Raphson algorithm to find the root of the
> function"f".
Can you Chris provide us an example with concrete functions? To us, it
would be easier to think about a concrete example.
Paul
__
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] numerical integration of a ftn of 2 variables
Chris Rhoads wrote: > > > I wish to find the root of a function of two variables that is defined by > an integral which must be > evaluated numerically. > > So the problem I want to solve is of the form: Find k such that f(k)=0, > where f(y) = int_a^b > g(x,y) dx. Again, the integral > involved must be done numerically. > > g <- function(x,p) 2*(x-p) Find value of p such that integral_0^1 g(x,p)dx is zero f <- function(p) integrate(g,0,1,p)$value uniroot(f, c(0,1)) $root should be 0.5 which is what is expected here. Berend Hasselman -- View this message in context: http://www.nabble.com/numerical-integration-of-a-ftn-of-2-variables-tp15578652p15618830.html Sent from the R help mailing list archive at Nabble.com. __ 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.
