[R] birthday problem (factorial limit)

Daniel Nordlund djnordlund at verizon.net
Sun Sep 28 21:40:36 CEST 2008


> -----Original Message-----
> From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org] On Behalf
> Of plummer at iarc.fr
> Sent: Sunday, September 28, 2008 11:35 AM
> To: ted.harding at manchester.ac.uk
> Cc: r-help at r-project.org; Uwe Ligges
> Subject: Re: [R] birthday problem (factorial limit)
> 
> Quoting "(Ted Harding)" <Ted.Harding at manchester.ac.uk>:
> 
> > On 28-Sep-08 17:51:55, Uwe Ligges wrote:
> > > Jörg Groß wrote:
> > >> Hi,
> > >> I tried to calculate the formula for the birthday problem
> > >> (the probability that at least two people out of a group of
> > >> n people share the same birthday)
> > >>
> > >> But the factorial-function allows me only to calculate
> > >> factorials up to 170.
> > >>
> > >> So is there a way to push that limit?
> > >>
> > >> to solve this formula:
> > >>
> > >> (factorial(365) / factorial((365-23))) / (365^23)
> > >
> > > Obviously you can easily rewrite this formula to:
> > >
> > > prod(343:365) / (365^23)
> > >
> > > or
> > >
> > > factorial(23) * choose(365, 23) / (365^23)
> > >
> > > Uwe Ligges
> > >
> > >> (n=23)
> >
> > I would put it in an even "safer" form:
> >
> > n <- 23
> > prod( ((365-(n-1)):365)/rep(365,n) )
> >
> > In other word: It evaluates
> >
> >  (343/365)*(344/365)* ... *(365/365)
> >
> > 365^N --> "Inf" if N > 120, whereas
> >
> >   n<-150
> >   prod( ((365-(n-1)):365)/rep(365,n) )
> > # [1] 2.451222e-16
> >
> > Best wishes,
> > Ted.
> 
> There is a pbirthday function and a qbirthday function in the
> stats package, which give approximate probability calculations
> for the birthday problem. The Examples section of the help page
> also contains Ted's more accurate solution.
> 
> Martyn

You can also use the lfactorial() or lgamma() functions as has been suggested in other discussions of factorials on this list.

exp((lfactorial(365) - lfactorial(365-23)) - 23*log(365))

or 
exp((lgamma(365 + 1) - lgamma(365-23 + 1)) - 23*log(365))


Hope this is helpful,

Dan

Daniel Nordlund
Bothell, WA USA 



More information about the R-help mailing list