[R] factorials

[Prof Brian Ripley]

On Tue, 15 Jan 2002, Damien Joly wrote:

> I'm a total newbie at using R, and so there probably
> is a better way to do this.  However, I couldn't find
> one, and so maybe this will help someone.
>
> I was calculating log-likelihoods using a multinomial
> model, and found that for large n, prod(n:1) wouldn't
> work to calculate factorials (e.g., prod(200:1) =
> Inf).  The below function calculates the natural log
> of a factorial (e.g. factorial(100) returns ln(100!)
> or in other words 100! = exp(factorial(100)).
>
> factorial<-function(p) {P<-array(c(p:1),dim=c(p,1))
> P<-log(P)
> return(sum(P))}

1) Factorials rapidly get large, so you need to work directly with logs.
In this case 200! is larger than the largest representable number (on a
IEC60559 machine, as all current R implementations seem to be).

2) log(n!) = lgamma(n+1), so use lgamma (in preference to sum(log(n:2)):
you did not need an array, especially not a 1-dimensional one, as R
objects are vectors).

3) R-help has discussed this fairly recently and you might like to refer
back to the archive.


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