[R] Product of 1 - probabilities
murdoch at stats.uwo.ca
Thu May 21 17:19:51 CEST 2009
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.
> 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.
> -or something else
> Any help greatly appreciated
> [[alternative HTML version deleted]]
> R-help at r-project.org mailing list
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
More information about the R-help