[R] Product of 1 - probabilities

[Duncan Murdoch]

On 5/21/2009 10:15 AM, Mark Bilton wrote:
> I am having a slight problem with probabilities.
> 
> To calculate the final probability of an event p(F), we can take the product of the chance that each independent event that makes p(F) will NOT occur.
> So...
> p(F) = 1- ( (1-p(A)) * (1-p(B)) * (1-p(C))...(1-p(x)) )
> 
> If the chance of an event within the product occurring remains the same, we can therefore raise this probability to a power of the number of times that event occurs.
> e.g. rolling a dice p(A) = 1/6 of getting a 1...
> p(F) = 1 - (1- (1/6))^z 
> p(F) = 1 - (1-p(A))^z tells us the probabiltity of rolling a 1 'at least once' in z number of rolls.
> 
> So then to R...
> 
> if p(A) = 0.01; z = 4; p(F) = 0.039
> 
> obviously p(F) > p(A)
> 
> however the problem arises when we use very small numbers e.g. p(B) = 1 * 10^-30
> R understands this value
> However when you have 1-p(B) you get something very close to 1 as you expect...but R seems to think it is 1.
> So when I calculate p(F) = 1 - (1-p(B))^z = 1 to the power anything equals 1 so p(F) = 0 and not just close to zero BUT zero.
> It doesn't matter therefore if z = 1*10^1000, the answer is still zero !!
> 
> Obviously therefore now p(F) < p(B)
> 
> Is there any solution to my problem, e.g.
> - is it a problem with the sum (-) ? ie could I change the number of bits the number understands (however it seems strange that it can hold it as a value close to 0 but not close to 1)
> -or should I maybe use a function to calculate the exact answer ?

The problem is that R maintains about 15-16 decimal places of accuracy, 
and to that accuracy, 1- 1.e-30 is equal to 1.  You need to find a 
different formula, e.g.

(1-x)^z = exp(z * log(1-x)) = exp(z * log1p(-x))

So if you have z = 1.e20, and p(A) = 1.e-30, then

1 - (1-1.e-30)^1.e20 can be calculated in R as

1 - exp(1.e20 * log1p(-1.e-30))

which works out to 1.e-10.

Duncan Murdoch

> -or something else
> 
> Any help greatly appreciated
> Mark
> -
> 
> 
> 
>   
> 
> _________________________________________________________________
> 
> 
> 	[[alternative HTML version deleted]]
> 
> ______________________________________________
> R-help at 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.