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