Hi there,
I have written a function to calculate factorials as follows:
fact <- function(x) {
recurse <- x > 1
x[!recurse] <- 1
if( any(recurse) ) {
y <- x[recurse]
x[recurse] <- y * fact( y - 1 )
}
x
}
I want to be able to do the famous birthday problem, which will involve
the computation of 365!, however it shall get cancelled out during the
computation. To make it more clear, I want to work out:
n <- 1:80
prob <- 1 - ( fact(365) / ( fact( 365 - n ) * 365 ^ n ) )
This, if works, should gives a vector of 80 probabilities, that given a
room of n people, 2 will share the same birthday.
But R gives me the following Error message:
Error in fact(y - 1) : evaluation is nested too deeply: infinite
recursion?
I assume that 365 is simply too large? Which should imply that my
factorial function is not written to the optimised way. Anyone can
offer a solution? Much appreciated.
Sincerely,
-------------------------------------------
Ko-Kang Wang
University of Auckland
Auckland 1005
New Zealand
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !) To: r-help-request at
stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
> Hi there, > > I want to be able to do the famous birthday problem, which will involve > the computation of 365!, however it shall get cancelled out during the > computation. To make it more clear, I want to work out: > n <- 1:80 > prob <- 1 - ( fact(365) / ( fact( 365 - n ) * 365 ^ n ) )Rule 1: If you find yourself calculating factorials you are probably doing it the wrong way. Rule 2: Same for determinants, but that's another issue. Sorry to spoil it for you, but I can't resist:> 1 - cumprod(365:285/365)[1] 0.000000000 0.002739726 0.008204166 0.016355912 0.027135574 [6] 0.040462484 0.056235703 0.074335292 0.094623834 0.116948178 [11] 0.141141378 0.167024789 0.194410275 0.223102512 0.252901320 [16] 0.283604005 0.315007665 0.346911418 0.379118526 0.411438384 [21] 0.443688335 0.475695308 0.507297234 0.538344258 0.568699704 [26] 0.598240820 0.626859282 0.654461472 0.680968537 0.706316243 [31] 0.730454634 0.753347528 0.774971854 0.795316865 0.814383239 [36] 0.832182106 0.848734008 0.864067821 0.878219664 0.891231810 [41] 0.903151611 0.914030472 0.923922856 0.932885369 0.940975899 [46] 0.948252843 0.954774403 0.960597973 0.965779609 0.970373580 [51] 0.974431993 0.978004509 0.981138113 0.983876963 0.986262289 [56] 0.988332355 0.990122459 0.991664979 0.992989448 0.994122661 [61] 0.995088799 0.995909575 0.996604387 0.997190479 0.997683107 [66] 0.998095705 0.998440043 0.998726391 0.998963666 0.999159576 [71] 0.999320753 0.999452881 0.999560806 0.999648644 0.999719878 [76] 0.999777437 0.999823779 0.999860955 0.999890668 0.999914332 [81] 0.999933109 Hmm. Better than even chances with 23 people. Amazing... :-) -- Bill Venables, Statistician, CMIS Environmetrics Project CSIRO Marine Labs, PO Box 120, Cleveland, Qld, AUSTRALIA. 4163 Tel: +61 7 3826 7251 Email: Bill.Venables at cmis.csiro.au Fax: +61 7 3826 7304 http://www.cmis.csiro.au/bill.venables/ -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html Send "info", "help", or "[un]subscribe" (in the "body", not the subject !) To: r-help-request at stat.math.ethz.ch _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
On Tue, 25 Apr 2000, Ko-Kang Wang wrote:> Hi there, > > I have written a function to calculate factorials as follows: > > fact <- function(x) { > recurse <- x > 1 > x[!recurse] <- 1 > if( any(recurse) ) { > y <- x[recurse] > x[recurse] <- y * fact( y - 1 ) > } > x > } > > > I want to be able to do the famous birthday problem, which will involve > the computation of 365!, however it shall get cancelled out during the > computation. To make it more clear, I want to work out: > n <- 1:80 > prob <- 1 - ( fact(365) / ( fact( 365 - n ) * 365 ^ n ) ) > > This, if works, should gives a vector of 80 probabilities, that given a > room of n people, 2 will share the same birthday. > > But R gives me the following Error message: > Error in fact(y - 1) : evaluation is nested too deeply: infinite > recursion? > > I assume that 365 is simply too large? Which should imply that my > factorial function is not written to the optimised way. Anyone can > offer a solution? Much appreciated.You _could_ increase the iteration limit with options(expressions=5000), say. But, you could also use the fact that n! = gamma(n+1) so> gamma(366)[1] Inf> lgamma(366)[1] 1792.332 and your answer would overflow the machine arithmetic, at least on 32-bit machines, as the answer is about 10^778. As Bill Venables points out, you want to do the cancellation algebraically. My reason for chipping is to note that you _never_ need to calculate a factorial recursively, and that calculations like this are often best done on log scale. You can do the birthday problem that way n <- 1:50 1 - exp(lgamma(366) - lgamma(366-n) - n*log(365)) agrees with Bill's answer up to 12dp or so. -- Brian D. Ripley, ripley at stats.ox.ac.uk Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272860 (secr) Oxford OX1 3TG, UK Fax: +44 1865 272595 -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html Send "info", "help", or "[un]subscribe" (in the "body", not the subject !) To: r-help-request at stat.math.ethz.ch _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._