[R] Factorials (was Recursive Computation in R)

[Prof Brian D Ripley]

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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._