

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]]

