Try the following:
 
sech<-function(X) 2/(exp(-X)+exp(X))
your.integrand <-  function(X) {0.5*cos(theta)*exp(X*theta)*sech(pi*X/2)}

 
my.integrand <- function(X) {0.5 * cos(theta) * 2 / (exp(-pi*X/2 - X*theta)
+ exp(pi*X/2 - X*theta)) }
 
theta <-  -1
my.integrand(-800)
your.integrand(-800)
 
Do you see the problem?
 
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@jhmi.edu

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

 

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

 

  _____  

From: xia wang [mailto:wxialuck@hotmail.com] 
Sent: Tuesday, August 26, 2008 12:59 PM
To: Ravi Varadhan
Cc: r-help@r-project.org
Subject: RE: [R] Problem with Integrate for NEF-HS distribution


Thanks so much. It works well in my MCMC sampling.  May I know why
re-writing the integrand can solve the problem? I have been thinking it was
the skewness of the distribution brought the error.  Thanks!

Xia



> From: RVaradhan@jhmi.edu
> To: wxialuck@hotmail.com; r-help@r-project.org
> Subject: RE: [R] Problem with Integrate for NEF-HS distribution
> Date: Tue, 26 Aug 2008 11:59:44 -0400
> 
> 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@jhmi.edu
> 
> Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html
> 
> 
> 
>
----------------------------------------------------------------------------
> --------
> 
> 
> -----Original Message-----
> From: r-help-bounces@r-project.org [mailto:r-help-bounces@r-project.org]
On
> Behalf Of xia wang
> Sent: Tuesday, August 26, 2008 1:01 AM
> To: r-help@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@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.


  _____  

See what people are saying about Windows Live. Check out featured posts.
Check It Out!
<http://www.windowslive.com/connect?ocid=TXT_TAGLM_WL_connect2_082008>  

	[[alternative HTML version deleted]]

