[R] Problem with Integrate for NEF-HS distribution

[Ravi Varadhan]

Hi,

Simply re-write your integrand as follows:


integrand <- function(X) {0.5 * cos(theta) * 2  / (exp(-pi*X/2 - X*theta) +
exp(pi*X/2 - X*theta)) }

Now the problem goes away.

> theta <- -1
> integrate(integrand, -Inf, 1)
0.9842315 with absolute error < 1.2e-05

This would work for all theta such that abs(theta) < -pi/2.

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 xia wang
Sent: Tuesday, August 26, 2008 1:01 AM
To: r-help at r-project.org
Subject: [R] Problem with Integrate for NEF-HS distribution



I need to calcuate the cumulative probability for the Natural Exponential
Family - Hyperbolic secant distribution with a parameter theta between -pi/2
and pi/2. The integration should be between 0 and 1 as it is a probability. 

The function "integrate" works fine when the absolute value of theta is not
too large.  That is, the NEF-HS distribution is not too skewed.  However,
once the theta gets large in absolute value, such as -1 as shown in the
example below, "integrate" keeps give me error message for "non-finite
function" when I put the lower bound as -Inf.  I suspect that it is caused
by the very heavy tail of the distribution. 

Is there any way that I can get around of this and let "integrate" work for
the skewed distribution?  Or is there any other function for integrating in
R-package? Thanks a lot for your advice in advance!


> theta<--1
> sech<-function(X) 2/(exp(-X)+exp(X))
> integrand <- function(X) {0.5*cos(theta)*exp(X*theta)*sech(pi*X/2)}

> integrate(integrand, -3,1)
0.8134389 with absolute error < 7e-09
> integrate(integrand, -10,1)
0.9810894 with absolute error < 5.9e-06
> integrate(integrand, -15,1)
0.9840505 with absolute error < 7e-09
> integrate(integrand, -50,1)
0.9842315 with absolute error < 4.4e-05
> integrate(integrand, -100,1)
0.9842315 with absolute error < 3.2e-05
> integrate(integrand, -Inf,1)
Error in integrate(integrand, -Inf, 1) : non-finite function value


Xia
_________________________________________________________________
Be the filmmaker you always wanted to be-learn how to burn a DVD with
WindowsR.

______________________________________________
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.