2
$\begingroup$

I want to calculate ${n \choose k}/2^n$ for moderate $n$ and $k$. In Matlab, use nchoosek(n,k) with $n=60$ and $k=30$ will give a warning: "Warning: Result may not be exact. Coefficient is greater than 9.007199e+15 and is only accurate to 15 digits ". I understand that $60\choose30$ is a huge number, but what I really need is ${60\choose30}/2^{60}$, which is around 0.1026.

My question is, is there a way to compute ${n \choose k}/2^n$ without sacrificing numerical precision in Matlab?

$\endgroup$
  • 1
    $\begingroup$ What exactly is the problem? Matlab still gives you a result that is as precise as expressible in 64bit floating point and that can be used in further computations. You will not get more accuracy in floating point precision by changing the exponent of 2 in that number by -60, i.e., by dividing by $2^{60}$. So is this really about how to suppress this warning or do you need the full accuracy of this number? Then you need some big integer or big decimal package. $\endgroup$ – Lutz Lehmann Mar 15 '14 at 18:11
5
$\begingroup$

If you have MATLAB's Symbolic Math Toolbox installed, then it is just a matter of writing:

evalin(symengine, 'binomial(60, 30) / 2^60')

Alternatively, you could write your own version of nchoosek (see this) using multiple precision arithmetic (available in MATLAB through third party toolboxes like this one). You may also consider writing your own code for addition of arbitrary precision integer arithmetic (this is not too hard).

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ How do I get floating point value out of this (I'm new to Matlab's symbolic math toolbox... I have been using Mathematica exclusively for symbolic works...) $\endgroup$ – wdg Mar 16 '14 at 17:13
  • $\begingroup$ Use double(evalin(...)) $\endgroup$ – Juan M. Bello-Rivas Mar 16 '14 at 19:23
4
$\begingroup$

I'd take the log of your expression, calculate the log of your expression using built-in functions that are well-behaved (i.e., don't underflow or overflow), and then exponentiate at the end.

\begin{align} {n \choose k} = \frac{n!}{k!(n-k)!} \end{align}

so

\begin{align} \log\left(\frac{{n \choose k}}{2^{n}}\right) = \log{n!} - \log{k!} - \log{(n - k)!} - n \log{2}. \end{align}

You can calculate $\log{n!}$ via gammaln(n+1), since $n!$ is $\Gamma(n+1)$, where $\Gamma$ is the Gamma function. (Don't calculate $\log{n!}$ as log(gamma(n+1)), since it won't avoid overflow issues!) This method should be better behaved. You could also get a quick estimate by using Stirling's approximation to $\log{n!}$, which is $n \log{n} - n$; I suspect you want something more accurate, but Stirling's approximation will give you a reasonable idea as to whether or not your exponential operation at the end will overflow.

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ This suffers from some mild cancellation. For $n=200$ and $k=100$, for instance, the last two digits are inexact. My approach below should work in full precision instead, as the naive nchoosek(200,100)/2^200, despite the warning (result compared with Wolfram Alpha). $\endgroup$ – Federico Poloni Mar 15 '14 at 18:13
  • $\begingroup$ As I understand the matlab philosophy, if you make no extra effort, all results are in 64bit double. Why should this method give a more relatively accurate double result than the original nchoosek function, or its logarithm, in view of the fact that the exact result exceeds the precision of the double format? $\endgroup$ – Lutz Lehmann Mar 15 '14 at 18:30
  • $\begingroup$ You're both right; this approach mainly avoids overflow at the expense of mild cancellation. $\endgroup$ – Geoff Oxberry Mar 15 '14 at 19:05
3
$\begingroup$

If utmost speed is not a concern, I'd go for rewriting nchoosek interspersing the divisions with the computation so that the temporary values stay bounded. What follows is a sample implementation.

function accumulator=dividedbinomial(n,k)
  num=n; %factor to multiply at the numerator
  den=k; %factor to multiply at the denominator
  powers=n; %powers of 2 still left to divide
  accumulator=1; 
  for i=1:k
    accumulator=accumulator*n/k;
    n=n-1;
    k=k-1;
    %divides now by enough powers of two so that the result stays below 1
    while accumulator>1 && powers>0 
      accumulator=accumulator/2;
      powers=powers-1;
    end
  end
  %divides by the remaining powers of two
  while powers>0
    accumulator=accumulator/2;
    powers=powers-1;
  end
end
| cite | improve this answer | |
$\endgroup$
  • 1
    $\begingroup$ There are 30, 15, 7, 3, 1 factors in $60!$ that are divisible by 2, 4, 8, 16, and 32 giving 56 as the power of 2 in 60! and 26 as the power of 2 in 30!. $\binom{60}{30}$ contains thus 2 with multiplicity 56-2*26=4, the remaining 56 divisions by 2 lead to a fractional part, exceeding the precision of 64bit floating point numbers by 3 bits, and that only if one forgets about all the other factors. In what way do you expect that your procedure, hopefully computed in floating point and not integer arithmetic, gives more relative accuracy than the result of the original nchoosek function? $\endgroup$ – Lutz Lehmann Mar 15 '14 at 18:24
  • $\begingroup$ @LutzL I don't; the two should be equivalent, apart from the fact that this version does not overflow (try n=2000,k=1000 with nchoosek, for instance). The result is not an integer, so you will need to restrict to double precision anyway. The relative error should be bounded by $(2k+n)*eps$, and in practice this seems to return all correct digits even for n=2000, k=1000. I don't think that there is a way to make it more accurate than this without resorting to higher-precision arithmetic. $\endgroup$ – Federico Poloni Mar 15 '14 at 18:30
  • $\begingroup$ Yes, if you expand the question in this direction, then your procedure ensures to get a result at all. The answer of Juan M. Bello Rivas shows how to invoke higher-precision arithmetic. $\endgroup$ – Lutz Lehmann Mar 15 '14 at 18:33
  • $\begingroup$ @LutzL Yes, I agree completely with your analysis: in floating point arithmetic clearly we cannot do better than what nchoosek already does (assuming it works correctly), especially so since dividing by 2 "comes for free" in that representation. We can only have marginal improvements such as no overflow. $\endgroup$ – Federico Poloni Mar 15 '14 at 19:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.