Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

R: How can I calculate large numbers in n-choose-k? [duplicate]

For a class assignment, I need to create a function that calculates n Choose k. I did just that, and it works fine with small numbers (e.g. 6 choose 2), but I'm supposed to get it work with 200 choose 50, where it naturally doesn't. The answer is too large and R outputs NaN or Inf, saying:

> q5(200, 50)
[1] "NaN"
Warning message:
In factorial(n) : value out of range in 'gammafn'

I tried using logs and exponents, but it doesn't cut it.

q5 <- function (n, k) {
  answer <- log(exp( factorial(n) / ( (factorial(k)) * (factorial(n - k)) )))
  paste0(answer)
}
like image 698
Ronny Efronny Avatar asked Nov 10 '16 11:11

Ronny Efronny


2 Answers

The answer to the actual question is that R cannot show numbers it cannot represent, and some of the terms in your equation are too big to represent. So it fails. However there are approximations to factorial that can be used - they work with logarithms which get big a lot slower.

The most famous one, Sterling's approximation, was not accurate enough, but the Ramanujan's approximation came to the rescue :)

ramanujan <- function(n){
  n*log(n) - n + log(n*(1 + 4*n*(1+2*n)))/6 + log(pi)/2
}

nchoosek <- function(n,k){
  factorial(n)/(factorial(k)*factorial(n-k))
} 

bignchoosek <- function(n,k){
  exp(ramanujan(n) - ramanujan(k) - ramanujan(n-k))
}

nchoosek(20,5)
# [1] 15504

bignchoosek(20,5)
# [1] 15504.06


bignchoosek(200,50)
# [1] 4.538584e+47
like image 192
Mike Wise Avatar answered Sep 28 '22 11:09

Mike Wise


You can try this too:

q5 <- function (n, k) {
  # nchoosek = (n-k+1)(n-k+2)...n / (1.2...k)
  return(prod(sapply(1:k, function(i)(n-k+i)/(i))))
}

q5(200, 50)
#[1] 4.538584e+47

or in log domain

q5 <- function (n, k) {
  # ln (nchoosek) = ln(n-k+1) + ln(n-k+2) + ...+ ln(n) - ln(1) -ln(2) - ...- ln(k)
  return(exp(sum(sapply(1:k, function(i)(log(n-k+i) - log(i))))))
}
q5(200, 50)
#[1] 4.538584e+47
like image 28
Sandipan Dey Avatar answered Sep 28 '22 10:09

Sandipan Dey