多项式系数的高效 Matlab 实现答案

【问题标题】：Efficient Matlab implementation of Multinomial Coefficient多项式系数的高效 Matlab 实现
【发布时间】：2014-03-10 13:11:55
【问题描述】：

我要计算多项式系数：

满足的地方n=n0+n1+n2

这个算子的Matlab实现可以很方便的在函数中完成：

function N = nchooseks(k1,k2,k3)
    N = factorial(k1+k2+k3)/(factorial(k1)*factorial(k2)*factorial(k3)); 
end

但是，当索引大于 170 时，阶乘将是无限的，在某些情况下会生成 NaN，例如180!/(175! 3! 2!) -> Inf/Inf-> NaN.

在其他帖子中，他们已经解决了 C 和 Python 的溢出问题。

在 C: 的情况下，你可以把所有的阶乘组成列表，然后找到列表中所有数字的素因式分解，然后将顶部的所有数字与底部的数字相消，直到数量完全减少”。
对于 Python：“利用阶乘 (n) = gamma(n+1) 的事实，使用 gamma 函数的对数并使用加法而不是乘法，减法而不是除法” 。

第一个解决方案似乎非常慢，所以我尝试了第二个选项：

function N = nchooseks(k1,k2,k3)
    N = 10^(log_gamma(k1+k2+k3)-(log_gamma(k1)+log_gamma(k2)+log_gamma(k3))); 
end
function y = log_gamma(x),  y = log10(gamma(x+1));  end

我将原始和 log_gamma 实现与以下代码进行比较：

% Calculate
N=100; err = zeros(N,N,N);
for n1=1:N,
    for n2=1:N,
        for n3=1:N,
            N1 = factorial(n1+n2+n3)/(factorial(n1)*factorial(n2)*factorial(n3)); 
            N2 = 10^(log10(gamma(n1+n2+n3+1))-(log10(gamma(n1+1))+log10(gamma(n2+1))+log10(gamma(n3+1)))); 
            err(n1,n2,n3) = abs(N1-N2); 
        end
    end
end
% Plot histogram of errors
err_ = err(~isnan(err));
[nelements,centers] = hist(err_(:),1e2);
figure; bar(centers,nelements./numel(err_(:)));

但是，某些情况下的结果略有不同，如下面的直方图所示。

因此，我应该假设我的实现是正确的，还是数字错误不能证明数字分歧是合理的？

【问题讨论】：

要计算伽马的日志，您可以直接使用gammaln
相关：stackoverflow.com/questions/15301885/…

标签： matlab overflow factorial multinomial

【解决方案1】：

为什么不使用这个？速度快，不会溢出：

N = prod([1:n]./[1:n0 1:n1 1:n2]);

【讨论】：

我已经尝试过您的巧妙解决方案，但我得到了类似的错误分布。因此，我假设这三个实现是相似的。但是，您的解决方案是最快的（factorial:57s, log_gamma:15, product:6s）。因此，我为您提供解决方案。谢谢你:)

【解决方案2】：

很抱歉恢复旧帖子，但对于未来的搜索者，您几乎肯定应该将多项式系数写为二项式系数的乘积，并使用内置方法来计算二项式系数（或者使用 Pascal 的方法编写您自己的三角形或其他方法）。相关公式出现在Wikipedia section on multinomial coefficients 的第一段。（我会在这里写，但似乎没有渲染 LaTeX 的方法。）

这种方法的另一个好处是它尽可能好地解决溢出问题，因为因子都是整数。计算多项式系数时，不需要除法。

【讨论】：

【解决方案3】：

使用@jemidiah 提供的提示，

这是代码

function c = multicoeff (k), 
    c = 1; 
    for i=1:length(k), 
      c = c* bincoeff(sum(k(1:i)),k(i)); 
    end; 
end

以及一些使用示例：

octave:88> multicoeff([2 2 2])
ans =  90
octave:89> factorial(6)/(factorial(2)*factorial(2)*factorial(2))
ans =  90
octave:90> multicoeff([5 4 3])
ans =  27720
octave:91> factorial(12)/(factorial(5)*factorial(4)*factorial(3))
ans =  27720

【讨论】：

【解决方案4】：

另一种方法是使用 Yannis Manolopoulos 迭代法。假设我们有一个带有多项式条目的向量k。

function N = multicoeff (k),
  n=sum(k); 
  [_,imax]=max(k); 
  num=[n:-1:n-k(imax)-1]; 
  den=[]; k(imax)=[]; 
  for i=1:length(k), den=[den 1:k(i)]; endfor; 
  N=prod(num./den);
endfunction

例子

octave:2> k = [5 4 3];
octave:3> multicoeff (k)
ans =  27720

参考：雅尼斯·马诺洛普洛斯。二项式系数计算。 ACM SIGCSE 公告，34(4):65，2002 年 12 月。doi：10.1145/820127.820168。网址https: //doi.org/10.1145/820127.820168。

【讨论】：