3

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

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

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

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

4

3 回答 3

3

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

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);
于 2013-11-21T00:04:22.533 回答
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对数更糟糕吗？或相反亦然？

于 2013-11-21T00:56:35.910 回答
0

这仅适用于 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
于 2013-11-20T23:51:05.647 回答