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 <[EMAIL PROTECTED]> wrote:
> From: xia wang <[EMAIL PROTECTED]>
> Subject: RE: [R] Problem with Integrate for NEF-HS distribution
> To: [EMAIL PROTECTED], [EMAIL PROTECTED]
> 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: [EMAIL PROTECTED]
> > Subject: Re: [R] Problem with Integrate for NEF-HS
> distribution
> > To: [EMAIL PROTECTED]
> >
> > Shouldn't you define
> > integrand <- function(X,theta) ......
> > and not
> > integrand <- function(X) .....
> >
> >
> > --- On Tue, 26/8/08, xia wang
> <[EMAIL PROTECTED]> wrote:
> >
> > > From: xia wang <[EMAIL PROTECTED]>
> > > 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®.
> http://clk.atdmt.com/MRT/go/108588797/direct/01/
______________________________________________
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.
