Got it. Thanks so much!
Xia
From: [EMAIL PROTECTED]
To: [EMAIL PROTECTED]
CC: r-help@r-project.org
Subject: RE: [R] Problem with Integrate for NEF-HS distribution
Date: Tue, 26 Aug 2008 14:50:18 -0400
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: [EMAIL PROTECTED]
Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html
------------------------------------------------------------------------------------
From: xia wang [mailto:[EMAIL PROTECTED]
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: [EMAIL PROTECTED]
> To:
[EMAIL PROTECTED]; 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: [EMAIL PROTECTED]
>
> Webpage:
>
>
>
>
----------------------------------------------------------------------------
>
--------
>
>
> -----Original Message-----
> From:
[EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] 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.
eck It Out!
_________________________________________________________________
[[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.
