[R] Types of quadrature

Thu Mar 13 17:29:17 CET 2008

David,

The problem is with 1 - pghyp(.).  Here is a better way to compute your
omega - I first compute a "complementary" pghyp, which is 1 - pghyp, and
then use this to compute the numerator.  The denominator is okay as it is.

pghyp.c <- function(x) sapply(x, function(x){integrate(function(x)dghyp(x),
lower=x, upper=Inf, subdivisions=1000, rel.tol=1.e-07)$val})

L <- 1

int.num <- integrate(function(x)pghyp.c(x), lower=L, upper=Inf)$val

int.denom <- integrate(pghyp, lower = -Inf, upper = L)$val

> int.num
[1] 0.08397642

> int.denom
[1] 1.083976

You should contact the package developer to provide an option for computing
1 - pghyp by specifying an option such as lower.tail=FALSE, as it is done
with pnorm().  This would also solve your problem.

Hope this helps,
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: rvaradhan at jhmi.edu

Webpage:  http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html

----------------------------------------------------------------------------
--------

-----Original Message-----
From: Lüthi David (luda) [mailto:luda at zhaw.ch] 
Sent: Thursday, March 13, 2008 12:16 PM
To: Ravi Varadhan
Subject: AW: [R] Types of quadrature

Hi Ravi

Thanks again. The integral exists because the ghyp density decays
exponentially for x -> \infty. The problem is that pghyp (which performs
numerical integration of the density dghyp) gives 0 for large x (e.g. 20000)
which results in (1 - pghyp) = 1 for large x and this of course diverges. I
know that there is a package rjacobi and also cwhmisc with function
adaptsim, however, both of them seem to be inappropriate.
Adaptsim does not accept infinite boundaries and rjacobi does need
quadrature nodes...

Best regards,
David

-----Ursprüngliche Nachricht-----
Von: Ravi Varadhan [mailto:rvaradhan at jhmi.edu] 
Gesendet: Donnerstag, 13. März 2008 16:53
An: Lüthi David (luda)
Betreff: RE: [R] Types of quadrature

Hi David,

As integrate() says, your integrals for ghyp are probably divergent.  Do you
know any theoretical result that says that the ratio you are computing for
ghyp actually exists?  

Here is how you would integrate in your pnorm() example:

> integrate(function(x) pnorm(x, lower.tail=FALSE), lower=1, upper=Inf)

0.08331547 with absolute error < 1.5e-07

>

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: rvaradhan at jhmi.edu

Webpage:  http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html

----------------------------------------------------------------------------
--------

-----Original Message-----
From: Lüthi David (luda) [mailto:luda at zhaw.ch] 
Sent: Thursday, March 13, 2008 4:41 AM
To: Ravi Varadhan
Subject: AW: [R] Types of quadrature

Dear Prof. Ripley, Dear Prof. Varadhan

Forgive me for being unprecise. Let me demonstrate my problem in an example:

Assume the cumulative distribution function has to be integrated
numerically:

scalar.pnorm <- function(x){
  return(integrate(dnorm, -Inf, x)$value)
}

.pnorm <- Vectorize(scalar.pnorm)

integrate(function(x)1 - .pnorm(x), lower = 1, upper = Inf)$value

The final usage would be:

library(ghyp)

omega.ghyp <- function(L, object = ghyp(), ...){
  int.num <- integrate(function(x, object, ...)1 - pghyp(x, object, ...),
object = object, lower = L, upper = Inf, ...) 
  int.denom <- integrate(pghyp, object = object, lower = -Inf, upper = L,
...)
  return(int.num$value / int.denom$value)
}

object <- ghyp()
omega.ghyp(1, object)

I do not see a way to avoid this or reformulate those expressions!?

Many thanks and best regards,

David

-----Ursprüngliche Nachricht-----
Von: Ravi Varadhan [mailto:rvaradhan at jhmi.edu] 
Gesendet: Mittwoch, 12. März 2008 19:21
An: Lüthi David (luda); r-help at r-project.org
Betreff: RE: [R] Types of quadrature

Hi,

Why do you need an extension of integrate()?  integrate() is adaptive - it
uses an adaptive Gauss-Kronrod quadrature.

You can specify Inf and -Inf as upper and lower limits, resp., in
integrate().  In fact, this is what the help page recommends, and it also
discourages the use of a large number as a surrogate for Inf.

What is the specific problem or distribution that you are having trouble
with in using integrate()?

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: rvaradhan at jhmi.edu

Webpage:  http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html

----------------------------------------------------------------------------
--------

-----Original Message-----
From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org] On
Behalf Of Lüthi David (luda)
Sent: Wednesday, March 12, 2008 1:25 PM
To: r-help at r-project.org
Subject: [R] Types of quadrature

Dear R-users

I would like to integrate something like \int_k^\infty (1 - F(x)) dx, where
(.) is a cumulative distribution function. As mentioned in the "integrate"
help-page: integrate(dnorm,0,20000) ## fails on many systems. This does not
happen for an adaptive Simpson or Lobatto quadrature (cf. Matlab). Even
though I am hardly familiar with numerical integration the implementation
seems to be fairly straightforward. 

My questions: 
- Is this extension of the function "integrate" planned for upcoming
versions of R? 
- Do there exist packages / workarounds?

I'm using R 2.6.2 on Windows and the reason why I want to integrate such an
expression is for the sake to compute the performance measure "Omega" for
financial securities.

Best regards,
David

--
David Lüthi
idp - Institute of Data Analysis and Process Design
Zurich University of Applied Sciences
Postfach 805
CH-8401 Winterthur

E-mail: david.luethi at zhaw.ch
Phone: 058 934 78 03
http://www.idp.zhaw.ch 
--

______________________________________________
R-help at 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.