scipy.misc.comb, returning n choose k, is implemented using the gammaln function. Is there a function that stays in log space? I see there is no scipy.misc.combln or any similar. It is trivial to implement myself, but it would be convenient if it were already in a package somewhere. I don't see it in scipy.misc, and it just feels wasteful to convert to normal space and then back to log.
2 Answers
It's possible to use gammaln, but precision is lost in subtractions when N >> k.
This can be avoided via the relation to the beta function:
from numpy import log
from scipy.special import betaln
def binomln(n, k):
# Assumes binom(n, k) >= 0
return -betaln(1 + n - k, 1 + k) - log(n + 1)
Looking at the source code, it seems like you're right that it would be trivial to implement, but that it likely isn't implemented elsewhere in scipy. On the plus side, there's some error checking going on, so you could eliminate some if you do those checks elsewhere (this is similar to getting rid of the exponential). If you know you'll always give 0 <= k <= N, and each of k, N as an array, then it is down to:
from scipy import special
def chooseln(N, k)
return special.gammaln(N+1) - special.gammaln(N-k+1) - special.gammaln(k+1)
