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 Windows?.
Thanks. I revised the code but it doesn't make difference. The cause of the
error is just as stated in the ?integrate " If the function is
approximately constant (in particular, zero) over nearly all
its range it is possible that the result and error estimate may be seriously
wrong. " I have searched R-archive. It may not be feasible to solve the
problem by "rectangle approximation" as some postings suggested
because the integration is in fact within MCMC samplings as follows. The lower
bound is always -Inf. The upper bound and the value of theta changes for each
sampling in MCMC.
I tried to multiple the integrand by a large number ( like 10^16) and changes
the rel.tol. It does help for some range of theta but not all.
Xia
like.fun <- function(beta,theta) {
integrand <- function(X,theta) {0.5*cos(theta)*exp(X*theta)*sech(pi*X/2)}
upper<-x%*%beta
prob<-matrix(NA, nrow(covariate),1)
for (i in 1:nrow(covariate))
{prob[i]<-integrate(integrand,lower=-Inf,
upper=upper[i],theta=theta)$value
}
likelihood <- sum(log(dbinom(y,n,prob)))
return(likelihood)
}
> Date: Tue, 26 Aug 2008 00:49:02 -0700
> From: m_olshansky@yahoo.com
> Subject: Re: [R] Problem with Integrate for NEF-HS distribution
> To: wxialuck@hotmail.com
>
> Shouldn't you define
> integrand <- function(X,theta) ......
> and not
> integrand <- function(X) .....
>
>
> --- On Tue, 26/8/08, xia wang <wxialuck@hotmail.com> wrote:
>
> > From: xia wang <wxialuck@hotmail.com>
> > Subject: [R] Problem with Integrate for NEF-HS distribution
> > To: r-help@r-project.org
> > Received: Tuesday, 26 August, 2008, 3:00 PM
> > 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 Windows®.
> >
> > ______________________________________________
> > 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.
_________________________________________________________________
Be the filmmaker you always wanted to be—learn how to burn a DVD with Windows®.
[[alternative HTML version deleted]]
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.