Эффективная реализация Matlab полиномиального коэффициента

Я хочу рассчитать полиномиальный коэффициент:

введите здесь описание изображения

где это удовлетворено 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: "используйте тот факт, что factorial(n) = gamma(n+1), используйте логарифм гамма-функции и используйте сложения вместо умножения, вычитания вместо деления".

Первое решение кажется очень медленным, поэтому я попробовал второй вариант:

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_(:)));

Однако в некоторых случаях результаты немного отличаются, как показано на следующей гистограмме.

введите здесь описание изображения

Таким образом, должен ли я предполагать, что моя реализация верна, или числовая ошибка не оправдывает расхождение чисел?


matlab overflow factorial multinomial
person tashuhka    schedule 10.03.2014    source источник
comment
Чтобы вычислить журнал гаммы, вы можете использовать gammaln напрямую.   -  person Luis Mendo    schedule 10.03.2014
comment
Связано: stackoverflow.com/ вопросы/15301885/   -  person tashuhka    schedule 30.03.2014

Ответы (4)


arrow_upward
3
arrow_downward

Почему бы не использовать это? Это быстро и не страдает от переполнения:

N = prod([1:n]./[1:n0 1:n1 1:n2]);
person Luis Mendo    schedule 10.03.2014
comment
Я попробовал ваше отличное решение, и я получаю аналогичное распределение ошибок. Таким образом, я предполагаю, что три реализации похожи. Однако ваше решение является самым быстрым (factorial: 57s, log_gamma: 15, product: 6s). Таким образом, я принимаю решение. Спасибо :) - person tashuhka; 10.03.2014

arrow_upward
2
arrow_downward

Извините, что воскрешаю старый пост, но для будущих искателей вы почти наверняка должны просто написать свой полиномиальный коэффициент как произведение биномиальных коэффициентов и использовать встроенный метод для вычисления биномиальных коэффициентов (или написать свой собственный, используя треугольник Паскаля или другой метод). Соответствующая формула приведена в первом абзаце раздела Википедии о мультиномиальных коэффициентах. (Я бы написал это здесь, но, похоже, нет способа отобразить LaTeX.)

Еще одним преимуществом этого подхода является то, что он настолько хорош, насколько это возможно, с переполнением, поскольку все коэффициенты являются целыми числами. Нет внутренней необходимости делить при вычислении полиномиальных коэффициентов.

person Joshua P. Swanson    schedule 28.09.2015

arrow_upward
0
arrow_downward

используя совет, предоставленный @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
person LEo    schedule 08.04.2019

arrow_upward
0
arrow_downward

Другой подход заключается в использовании итеративного метода Янниса Манолопулоса. Предположим, у нас есть вектор 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 г. doi: 10.1145/820127.820168. URL-адрес http://https:%20//doi.org/10.1145/820127.820168.

person LEo    schedule 28.02.2020