[R] Factorials (was Recursive Computation in R)
Prof Brian D Ripley
ripley at stats.ox.ac.uk
Tue Apr 25 08:49:41 CEST 2000
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
More information about the R-help
mailing list