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: [EMAIL PROTECTED]
Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html
----------------------------------------------------------------------------
--------
-----Original Message-----
From: Lüthi David (luda) [mailto:[EMAIL PROTECTED]
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:[EMAIL PROTECTED]
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: [EMAIL PROTECTED]
Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html
----------------------------------------------------------------------------
--------
-----Original Message-----
From: Lüthi David (luda) [mailto:[EMAIL PROTECTED]
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:[EMAIL PROTECTED]
Gesendet: Mittwoch, 12. März 2008 19:21
An: Lüthi David (luda); r-help@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: [EMAIL PROTECTED]
Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html
----------------------------------------------------------------------------
--------
-----Original Message-----
From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On
Behalf Of Lüthi David (luda)
Sent: Wednesday, March 12, 2008 1:25 PM
To: r-help@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: [EMAIL PROTECTED]
Phone: 058 934 78 03
http://www.idp.zhaw.ch
--
______________________________________________
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.
