Get binomial coefficients

Question

In an attempt to vectorize a particular piece of Matlab code, I could not find a straightforward function to generate a list of the binomial coefficients. The best I could find was nchoosek, but for some inexplicable reason this function only accepts integers (not vectors of integers). My current solution looks like this:

mybinom = @(n) arrayfun(@nchoosek, n*ones(1,n), 1:n)

This generates the set of binomial coefficients for a given value of n. However, since the binomial coefficients are always symmetric, I know that I am doing twice as much work as necessary. I'm sure that I could create a solution that exploits the symmetry, but I'm sure that it would be at the expense of readability.

Is there a more elegant solution than this, perhaps using a Matlab function that I am not aware of? Note that I am not interested in using the symbolic toolbox.

Answer 1

If you want to minimize operations you can go along these lines:

n = 6;

k = 1:n;
result = [1 cumprod((n-k+1)./k)]

>> result

result =

     1     6    15    20    15     6     1

This requires very few operations per coefficient, because each cofficient is obtained exploiting the previously computed one.

You can reduce the number of operations by approximately half if you take into account the symmetry:

m1 = floor(n/2);
m2 = ceil(n/2);
k = 1:m2;
result = [1 cumprod((n-k+1)./k)];
result(n+1:-1:m1+2) = result(1:m2);

Answer 2

What about a modified version of Luis Mendo's solution - but in logarithms:

n = 1e4;
m1 = floor(n/2);
m2 = ceil(n/2);
k = 1:m2;

% Attempt to compute real value
out0 = [1 cumprod((n-k+1)./k)];
out0(n+1:-1:m1+2) = out0(1:m2);

% In logarithms
out1 = [0 cumsum((log(n-k+1)) - log(k))];
out1(n+1:-1:m1+2) = out1(1:m2);

plot(log(out0) - out1, 'o-')

The advantage of working with logarithms is that you can set n = 1e4; and still obtain a good approximation of the real value (nchoosek(1e4, 5e3) returns Inf and this is not a good approximation at all!).

EDIT following horchler's comment

You can use the gammaln function to obtain the same result but it's not faster. The two approximations seem to be quite different:

n = 1e7;
m1 = floor(n/2);
m2 = ceil(n/2);
k = 1:m2;

% In logarithms
tic
out1 = [0 cumsum((log(n-k+1)) - log(k))];
out1(n+1:-1:m1+2) = out1(1:m2);
toc
% Elapsed time is 0.912649 seconds.

tic
k = 0:m2;
out2 = gammaln(n + 1) - gammaln(k + 1) - gammaln(n - k + 1);
out2(n+1:-1:m1+2) = out2(1:m2);
toc
% Elapsed time is 1.020188 seconds.

tmp = out2 - out1;
plot(tmp, '.')

prctile(tmp, [0 2.5 25 50 75 97.5 100])
%   1.0e-006 *
%    -0.2217   -0.1462   -0.0373    0.0363    0.1225    0.2943    0.3846

Is adding three gammaln worse than adding n logarithms? Or viceversa?

Answer 3

This works for Octave only

You can use bincoeff function.
Example: bincoeff(5, 0:5)

EDIT : Only improvement I can think of goes like this. Maybe you already thought this trivial solution and didn't like it.

# Calculate only the first half
mybinomhalf = @(n) arrayfun(@nchoosek, n*ones(1,n/2+1), 0:n/2)

# pad your array symmetrically
mybinom = @(n) padarray(mybinomhalf(n), [0 n/2], 'symmetric', 'post')

# I couldn't test it and this line may not work