[R] Problem with Integrate for NEF-HS distribution

[Moshe Olshansky]

If you look at your sech(pi*x/2) you can write it as 
sech(pi*x/2) = 2*exp(pi*x/2)/(1 + exp(pi*x))
For x < -15, exp(pi*x) < 10^-20, so for this interval you can replace sech(pi*x/2) by 2*exp(pi*x/2) and so the integral from -Inf to -15 (or even -10 - depends on your accuracy requirements) can be computed analytically, so you are left with integral from -10 (or -15) to your upper bound and this should be all right for numerical integration.


--- On Wed, 27/8/08, xia wang <wxialuck at hotmail.com> wrote:

> From: xia wang <wxialuck at hotmail.com>
> Subject: RE: [R] Problem with Integrate for NEF-HS distribution
> To: m_olshansky at yahoo.com, r-help at lists.r-project.org
> Received: Wednesday, 27 August, 2008, 12:26 AM
> 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 at yahoo.com
> > Subject: Re: [R] Problem with Integrate for NEF-HS
> distribution
> > To: wxialuck at hotmail.com
> > 
> > Shouldn't you define
> > integrand <- function(X,theta) ......
> > and not
> > integrand <- function(X) .....
> > 
> > 
> > --- On Tue, 26/8/08, xia wang
> <wxialuck at hotmail.com> wrote:
> > 
> > > From: xia wang <wxialuck at hotmail.com>
> > > Subject: [R] Problem with Integrate for NEF-HS
> distribution
> > > To: r-help at 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 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.
> 
> _________________________________________________________________
> Be the filmmaker you always wanted to be—learn how to
> burn a DVD with Windows®.
> http://clk.atdmt.com/MRT/go/108588797/direct/01/