如何更有效地计算二项式的总和？

Question

我必须计算一个等式如下：

其中 k1,k2 给出。我正在使用 MATLAB 计算 P。我想我对上面的等式有一个正确的实现。但是，我的实施速度很慢。我认为问题出在二项式系数上。根据等式，我可以有一种有效的方法来加快时间吗？谢谢大家。

对于k1=150; k2=150; D=200;，需要11.6 秒

function main
warning ('off');
  function test_binom()
      k1=150; k2=150; D=200; P=0;
      for i=0:D-1
          for j=0:i
              if (i-j>k2||j>k1) 
                  continue;
              end
              P=P+nchoosek(k1,j)*nchoosek(k2,i-j)/nchoosek((k1+k2),i);          
          end 
      end
  end
f = @()test_binom(); 
timeit(f)
end

更新：对于测量时间，我发现nchoosek是计算时间长的原因。因此，我将函数重写如下

function re=choose(n, k)
    if (k == 0)
        re=1;
    else
        re=(n * choose(n - 1, k - 1)) / k;
    end
end

现在，计算时间减少了0.25秒。请问有什么更好的方法吗？

Answer 1

您可以将 nchoosek 的结果保存到 table 以防止重复计算函数，还提供了二项式系数的实现：

%binomial coefficients
function nk=nchoosek2(n, k)
    if n-k > k
        nk = prod((k+1:n) .* prod((1:n-k).^ (-1/(n-k))));
    else
        nk = prod((n-k+1:n) .* prod((1:k).^ (-1/k)) ) ;
    end
end
%function to store and retrieve results of nchoosek to/from a table
function ret = choose (n,k, D, K1, K2)
    persistent binTable = zeros(max([D+1,K1+K2+1]) , D+1);
    if binTable(n+1,k+1) == 0
        binTable(n+1,k+1) = nchoosek2(n,k);
    end
    ret = binTable(n+1,k+1);
end

function P = tst()
    P=0;k1=150; k2=150; D=200; P=0;
    choose(1,0,D,k1,k2);
    for i = 0:D-1
        for j = j=max(i - k2 , 0):min (i,k1-1)
            P=P+choose(k1,j)*choose(k2,i-j)/choose((k1+k2),i);
        end
    end
end

你的代码与 nchoosek2 相比： online demo

Answer 2

您可以对所有过程进行矢量化，这将使它变得超快而无需 mex。

首先是 nchoosek 函数：

function C = nCk(n,k)
% use smaller k if available
k(k>n/2) = n-k(k>n/2);
k = k(:);
kmat = ones(numel(k),1)*(1:max(n-k));
kmat = kmat.*bsxfun(@le,kmat,(n-k));
pw = bsxfun(@power,kmat,-1./(n-k));
pw(kmat==0) = 1;
kleft = ones(numel(k),1)*(min(k):n);
kleft = kleft.*bsxfun(@gt,kleft,k);
t = bsxfun(@times,kleft,prod(pw,2));
t (kleft==0) = 1;
C = prod(t,2);
end

然后beta和P计算：

function P = binomial_coefficient(k1,k2,D)
warning ('off','MATLAB:nchoosek:LargeCoefficient');
i_ind = nonzeros(triu(ones(D,1)*(1:D)))-1;
j_ind = nonzeros(tril(ones(D,1)*(1:D+1)).')-1;
valid = ~(i_ind-j_ind>=k2 | j_ind>=k1);
i_ind = i_ind(valid);
j_ind = j_ind(valid);
beta = @(ii,jj) nCk(k1,jj).*nCk(k2,ii-jj)./nCk((k1+k2),ii);
b = beta(i_ind,j_ind);
P = sum(b(:));
end

并且执行时间从 10.674 下降到 0.49696 秒。

编辑：

采用@rahnema1 的想法，我设法使它变得更快，对所有独特的 nCk 计算使用 table，因此其中 none 将完成更多一次。使用与上面相同的 nCk 函数，这就是新的 binomial_coefficient 函数的外观：

function P = binomial_coefficient(k1,k2,D)
warning ('off','MATLAB:nchoosek:LargeCoefficient');
i_ind = nonzeros(triu(ones(D,1)*(1:D)))-1;
j_ind = nonzeros(tril(ones(D,1)*(1:D+1)).')-1;
valid = ~(i_ind-j_ind>=k2 | j_ind>=k1);
i_ind = i_ind(valid);
j_ind = j_ind(valid);
ni = numel(i_ind);
all_n = repelem([k1; k2; k1+k2],ni); % all n's to be used in thier order
all_k = [j_ind; i_ind-j_ind; i_ind]; % all k's to be used in thier order
% get all unique sets of 'n' and 'k':
sets_tbl = unique([all_n all_k],'rows');
uq_n = unique(sets_tbl(:,1));
nCk_tbl = zeros([max(all_n) max(all_k)+1]);
% compute all the needed values of nCk:
for s = 1:numel(uq_n)
    curret_n = uq_n(s);
    curret_k = sets_tbl(sets_tbl(:,1)==curret_n,2);
    nCk_tbl(curret_n,curret_k+1) = nCk(curret_n,curret_k).';
end
beta = @(ii,jj) nCk_tbl(k1,jj+1).*nCk_tbl(k2,ii-jj+1)./nCk_tbl((k1+k2),ii+1);
b = beta(i_ind,j_ind);
P = sum(b(:));
end

现在，当运行仅需 0.01212 秒时，它不仅是超快代码，而且是 flying-talking-super-fast 代码！

Answer 3

答案有一个非常好的方法：创建一个 table 的值，您生成一次并稍后访问（以及其他一些巧妙的优化）。

我要改变的一件事是计算二项式系数的方式。如果您查看计算 nchoosek(n, k) 和 nchoosek(n, k+1) 的阶乘，您两次都在重新计算 n!，而对于 k+1，您正在重新计算 k! 和乘以 k+1。（与 (n-k)! 类似。）

我们可以根据 nchoosek(n, k).

的值迭代计算 nchoosek(n, k+1)，而不是每次都丢弃计算。

function L=combList(n, maxk)
% Create a vector of length maxk containing
%   [nchoosek(n, 1), nchoosek(n, 2), ..., nchoosek(n, maxk)]
% Note: nchoosek(n, 0) == nchoosek(n, n) == 1

assert(maxk<=n, 'maxk must be less than or equal to n');

L = zeros(1,maxk);
L(1) = n;                    % nchoosek(n, 1) == n
for k = 2:maxk
   L(k) = L(k-1)*(n-k+1)/k;
end

在您的程序中，您只需为 k1、k2 和 k1+k2 创建 3 个具有适当限制的列表，然后对这些列表进行索引以生成总和。

如何更有效地计算二项式的总和？

How to compute sum of binomial more efficiently?

performance

time

matlab

vectorization

mex