R help - Jun 2004 - Greater than 1 or less than 1?

I have problem evaluating the expression h = exp(x)/(exp(exp(x))-1) for 
large negative x. This expression is actually the probability
that y = 1 when y is a Poisson random variable truncated at 0, hence 
must satisfy 0 <= h <= 1. However, when
x < -18, I may get an h value that is  larger than 1 while the true 
value should be a quantity that is smaller than but very close to 1.
For example when x = -19, h = 1.00000000031174. I tried to use 
different form of h and none of
them give me an h value that is less than or equal to 1. I also tried 
to find patterns, but discovered none: for some x < -18, h is below 1; 
for others
h is greater than 1. Is there any trick that enables me to obtain 
theoretically correct results from this expression?

Thanks.

Zhen Chen, Ph.D.
Assistant Professor of Biostatistics
Department of Biostatistics & Epidemiology
University of Pennsylvania School of Medicine
423 Guardian Drive -- 625 Blockley Hall
Philadelphia, PA  19104-6021

voice:  215-573-8545
fax:    215-573-4865
email:  zchen at cceb.upenn.edu

"Zhen Chen" <zchen at cceb.upenn.edu> writes:
> I have problem evaluating the expression h = exp(x)/(exp(exp(x))-1)
> for large negative x. This expression is actually the probability
> that y = 1 when y is a Poisson random variable truncated at 0, hence
> must satisfy 0 <= h <= 1. However, when
> x < -18, I may get an h value that is  larger than 1 while the true
> value should be a quantity that is smaller than but very close to 1.
> For example when x = -19, h = 1.00000000031174. I tried to use
> different form of h and none of
> them give me an h value that is less than or equal to 1. I also tried
> to find patterns, but discovered none: for some x < -18, h is below 1;
> for others
> h is greater than 1. Is there any trick that enables me to obtain
> theoretically correct results from this expression?
expm1() should help...

What platform is this? I don't see it until x=-30:
> exp(-30)
[1] 9.357623e-14
> exp(exp(-30))-1
[1] 9.348078e-14

which is of course "wrong" since exp(x) > x + 1 mathematically.
However, notice that in double precision calculations the latter value
is only accurate up to 1e-15 or so. 

-- 
   O__  ---- Peter Dalgaard             Blegdamsvej 3  
  c/ /'_ --- Dept. of Biostatistics     2200 Cph. N   
 (*) \(*) -- University of Copenhagen   Denmark      Ph: (+45) 35327918
~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk)             FAX: (+45) 35327907

On Sun, 20 Jun 2004, Zhen Chen wrote:
> I have problem evaluating the expression h = exp(x)/(exp(exp(x))-1) for
> large negative x. This expression is actually the probability
> that y = 1 when y is a Poisson random variable truncated at 0, hence
> must satisfy 0 <= h <= 1. However, when

You would be better off using the functions provided for the Poisson
distribution.

Either
dpois(1,exp(x))/ppois(0.5,exp(x),lower.tail=FALSE)
or to get 1-h
ppois(1.5,exp(x),lower.tail=FALSE)/ppois(0.5,exp(x),lower.tail=FALSE)

	-thomas

> Date: Sun, 20 Jun 2004 18:32:46 -0400
> From: "Zhen Chen" <zchen at cceb.upenn.edu>
> 
> I have problem evaluating the expression h = exp(x)/(exp(exp(x))-1) for 
> large negative x. This expression is actually the probability
> that y = 1 when y is a Poisson random variable truncated at 0, hence 
> must satisfy 0 <= h <= 1. However, when
> x < -18, I may get an h value that is  larger than 1 while the true 
> value should be a quantity that is smaller than but very close to 1.
> For example when x = -19, h = 1.00000000031174. I tried to use 
> different form of h and none of
> them give me an h value that is less than or equal to 1. I also tried 
> to find patterns, but discovered none: for some x < -18, h is below 1; 
> for others
> h is greater than 1. Is there any trick that enables me to obtain 
> theoretically correct results from this expression?
> 
exp(x)/(exp(exp(x)) - 1) is less than 1 for me for x = -18 with R-1.9.0
on a Solaris 9 system, a NetBSD 2.0C system and a Windows XP system.

You could try exp(x)/expm1(exp(x)), but if the longer expression doesn't
work, then this may not either, though it does give slightly different
results on my systems.

E.g.:
> x = -19
> exp(x)/(exp(exp(x)) - 1) - 1
[1] -2.047911e-09
> exp(x)/(expm1(exp(x))) - 1
[1] -2.801398e-09

log(1 + x) and exp(x) - 1 are often treated specially in math libraries
for when x << 0.

Ray Brownrigg