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 belearn 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 belearn how to burn a DVD with Windows®.
[[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.
