From: Prof Brian Ripley (ripley@stats.ox.ac.uk)
Date: Wed 16 Jan 2002 - 17:45:27 EST
Message-id: <Pine.LNX.4.31.0201160634480.25451-100000@gannet.stats>
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@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@stat.math.ethz.ch _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
This archive was generated by hypermail 2.1.3 : Thu 17 Jan 2002 - 11:10:11 EST