获取二项式系数

Question

在尝试向量化一段特定的 Matlab 代码时，我找不到一个简单的函数来生成二项式系数的列表。我能找到的最好的是nchoosek，但由于某些莫名其妙的原因，这个函数只接受整数（不是整数向量）。我当前的解决方案如下所示：

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

这会为给定的n 值生成一组二项式系数。然而，由于二项式系数总是对称的，我知道我做的工作量是必要的两倍。我确信我可以创建一个利用对称性的解决方案，但我确信它会以牺牲可读性为代价。

有没有比这更优雅的解决方案，也许使用我不知道的 Matlab 函数？请注意，我对使用符号工具箱不感兴趣。

Answer 1

如果您想最小化操作，您可以按照以下方式进行：

n = 6;

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

>> result

result =

     1     6    15    20    15     6     1

这对每个系数只需要很少的操作，因为每个系数都是利用先前计算的系数获得的。

如果考虑到对称性，您可以将操作次数减少大约一半：

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

Luis Mendo 解决方案的修改版本怎么样 - 但以对数形式：

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-')

使用对数的优点是您可以设置 n = 1e4; 并且仍然可以获得实际值的良好近似值（nchoosek(1e4, 5e3) 返回 Inf 而这不是一个很好的近似值全部！）。

根据 horchler 的评论进行编辑

您可以使用gammaln 函数获得相同的结果，但速度并不快。这两个近似值似乎完全不同：

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

添加三个gammaln 是否比添加n 对数更糟糕？反之亦然？

Answer 3

这仅适用于 Octave

你可以使用 bincoeff 函数。
示例：bincoeff(5, 0:5)

编辑： 我能想到的唯一改进是这样的。也许你已经想到了这个微不足道的解决方案并且不喜欢它。

# 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