Re: [R] Problem with Integrate for NEF-HS distribution

2008-08-26 Thread 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 <[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
> > >
> __

Re: [R] Problem with Integrate for NEF-HS distribution

2008-08-26 Thread xia wang 

Got it.  Thanks so much!

Xia

From: [EMAIL PROTECTED]
To: [EMAIL PROTECTED]
CC: r-help@r-project.org
Subject: RE: [R] Problem with Integrate for NEF-HS distribution
Date: Tue, 26 Aug 2008 14:50:18 -0400










Try the following:
 
sech<-function(X) 
2/(exp(-X)+exp(X))
your.integrand <-  function(X) 
{0.5*cos(theta)*exp(X*theta)*sech(pi*X/2)}

 
my.integrand <- function(X) {0.5 * 
cos(theta) * 2 / (exp(-pi*X/2 - X*theta) + exp(pi*X/2 - X*theta)) 
}
 
theta <-  
-1
my.integrand(-800)
your.integrand(-800)
 
Do you see the problem?
 
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: [EMAIL PROTECTED]
Webpage:  http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html






From: xia wang [mailto:[EMAIL PROTECTED] 

Sent: Tuesday, August 26, 2008 12:59 PM
To: Ravi 
Varadhan
Cc: r-help@r-project.org
Subject: RE: [R] Problem 
with Integrate for NEF-HS distribution


Thanks so much. It works well in my MCMC sampling.  May I know 
why re-writing the integrand can solve the problem? I have been thinking it was 
the skewness of the distribution brought the error.  
Thanks!

Xia



> From: [EMAIL PROTECTED]
> To: 
[EMAIL PROTECTED]; r-help@r-project.org
> Subject: RE: [R] Problem with 
Integrate for NEF-HS distribution
> Date: Tue, 26 Aug 2008 11:59:44 
-0400
> 
> 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: [EMAIL PROTECTED]
> 
> Webpage: 

> 

> 
> 
> 

> 

> 
> 
> -Original Message-
> From: 
[EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On
> 
Behalf Of xia wang
> Sent: Tuesday, August 26, 2008 1:01 AM
> To: 
r-help@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@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-gui

Re: [R] Problem with Integrate for NEF-HS distribution

2008-08-26 Thread Ravi Varadhan 
Try the following:
 
sech<-function(X) 2/(exp(-X)+exp(X))
your.integrand <-  function(X) {0.5*cos(theta)*exp(X*theta)*sech(pi*X/2)}

 
my.integrand <- function(X) {0.5 * cos(theta) * 2 / (exp(-pi*X/2 - X*theta)
+ exp(pi*X/2 - X*theta)) }
 
theta <-  -1
my.integrand(-800)
your.integrand(-800)
 
Do you see the problem?
 
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: [EMAIL PROTECTED]

Webpage:  http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html

 




 

  _  

From: xia wang [mailto:[EMAIL PROTECTED] 
Sent: Tuesday, August 26, 2008 12:59 PM
To: Ravi Varadhan
Cc: r-help@r-project.org
Subject: RE: [R] Problem with Integrate for NEF-HS distribution


Thanks so much. It works well in my MCMC sampling.  May I know why
re-writing the integrand can solve the problem? I have been thinking it was
the skewness of the distribution brought the error.  Thanks!

Xia



> From: [EMAIL PROTECTED]
> To: [EMAIL PROTECTED]; r-help@r-project.org
> Subject: RE: [R] Problem with Integrate for NEF-HS distribution
> Date: Tue, 26 Aug 2008 11:59:44 -0400
> 
> 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: [EMAIL PROTECTED]
> 
> Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html
> 
> 
> 
>

> 
> 
> 
> -----Original Message-
> From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED]
On
> Behalf Of xia wang
> Sent: Tuesday, August 26, 2008 1:01 AM
> To: r-help@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@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.


  _  

See what people are saying about Windows Live. Check out featured posts.
Check It Out!
<http://www.windowslive.com/connect?ocid=TXT_TAGLM_WL_connect2_082008>  

[[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.

Re: [R] Problem with Integrate for NEF-HS distribution

2008-08-26 Thread xia wang 


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®.

[[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.

Re: [R] Problem with Integrate for NEF-HS distribution

2008-08-26 Thread xia wang 

Thanks so much. It works well in my MCMC sampling.  May I know why re-writing 
the integrand can solve the problem? I have been thinking it was the skewness 
of the distribution brought the error.  Thanks!

Xia

 

> From: [EMAIL PROTECTED]
> To: [EMAIL PROTECTED]; r-help@r-project.org
> Subject: RE: [R] Problem with Integrate for NEF-HS distribution
> Date: Tue, 26 Aug 2008 11:59:44 -0400
> 
> 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: [EMAIL PROTECTED]
> 
> Webpage:  http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html
> 
>  
> 
> 
> 
> 
> 
> -Original Message-
> From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On
> Behalf Of xia wang
> Sent: Tuesday, August 26, 2008 1:01 AM
> To: r-help@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@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.

_


[[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.

Re: [R] Problem with Integrate for NEF-HS distribution

2008-08-26 Thread Ravi Varadhan 
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: [EMAIL PROTECTED]

Webpage:  http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html

 





-Original Message-
From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On
Behalf Of xia wang
Sent: Tuesday, August 26, 2008 1:01 AM
To: r-help@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@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.

__
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.

[R] Problem with Integrate for NEF-HS distribution

2008-08-25 Thread xia wang 


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.