> -----Original Message----- > From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf > Of [EMAIL PROTECTED] > Sent: Sunday, September 28, 2008 11:35 AM > To: [EMAIL PROTECTED] > Cc: r-help@r-project.org; Uwe Ligges > Subject: Re: [R] birthday problem (factorial limit) > > Quoting "(Ted Harding)" <[EMAIL PROTECTED]>: > > > 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 ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
