[R] Product of 1 - probabilities

[(Ted Harding)]

On 21-May-09 14:45:20, Gabor Grothendieck wrote:
> There are several arbitrary precision packages available: gmp (an
> interface to the GNU multi-precision library on CRAN) and
> bc (an R interface to the bc arbitrary precision calculator):
> http://r-bc.googlecode.com
> 
> There are also packages providing R interfaces to two computer
> algebra systems and they both support not only arbitrary
> precision but also exact calculation:
> 
> http://rsympy.googlecode.com
> http://ryacas.googlecode.com
> 
>> library(bc)
>> 1 - (1-10^-75)^10
> [1] 0
>> bc("1 - (1-10^-75)^10")
> [1]
> ".0000000000000000000000000000000000000000000000000000000000000000000000
> 000100000000000000000000000000"

Thanks for pointing this out, Gabor! In particular, bc is something
I often use (but on its own). It is good to know that R can connect
with it.

However, while your call above will indeed return an adequate answer
for 1 - (1-10^-75)^10, and can be assigned numerically by, say,

  w <- as.numeric(bc("1 - (1-10^-75)^10"))

this still does not let Mark off the hook of verifying that he will,
in the end, get a sensible answer. If this calculation were part of
a more elaborate computation in which, say, some formula called for
(1-w), then that would still evaluate to 1, with possible dire results.

While plain algebraic expressions can easily be evaluated in bc, its
repertoire of mathematical functions is limited. No doubt the other
mathematical packages you mention have much richer resources; but for
anything complicated requiring their use I would probably go for
working externally with such a package. However, that may be a naive
view, since I have not tried the R interfaces!

Nevertheless, I would still strongly recommend doing some analysis
of the expressions in a computation when limits of precision are an
issue, since often a simple expression can be obtained which, in R,
would give adequate results.

Ted.


> On Thu, May 21, 2009 at 10:15 AM, Mark Bilton <mcbilton at hotmail.com>
> 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 ?
>> -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.
>>
> 
> ______________________________________________
> 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.

--------------------------------------------------------------------
E-Mail: (Ted Harding) <Ted.Harding at manchester.ac.uk>
Fax-to-email: +44 (0)870 094 0861
Date: 21-May-09                                       Time: 16:39:16
------------------------------ XFMail ------------------------------