Предыстория здесь
topic91034.html.
Для понимания вставлю решение с постоянной

. Для биномиального распределения с заданными

и

функция вероятности находится следующим образом:
n = 6;
p = .2;
pmf = zeros(1,n+1);
for i=0:n
pmf(i+1) = nchoosek(n,i)*p^i*(1-p)^(n-i);
end
Теперь просто интегрируем функцию вероятности в требуемом диапазоне. Например, чтобы найти

:
n=2;
disp( sum(pmf(n+1:end)) );
Векторизация в этом случае не составляет труда. Если необходимо найти функции распределения сразу для

серий экспериментов:
n = 6;
m = 600;
p = rand(1,m);
pmf = zeros(n+1,m);
for i=0:n
pmf(i+1,:) = nchoosek(n,i)*p.^i.*(1-p).^(n-i);
end
Вычисление биномиального коэффициента при помощи функции
nchoosek производится один раз за цикл, никаких проблем.
В случае переменной вероятности

необходимо для каждого

найти все возможные комбинации вероятностей

и

. На помощь приходит все та же функция
nchoosek, а точнее второй вариант ее использования, когда первый аргумент — вектор:
n = 6;
p = rand(1,n);
p_ = 1-p;
pmf = zeros(1,n+1);
for i=0:n
P = nchoosek(p, i);
P_ = flipud( nchoosek(p_,n-i) );
pmf(i+1) = sum( prod([P P_], 2) );
end
Первая и вторая строка цикла выдают все возможные комбинации

и

в количестве

и

раз соответственно. Хотелось бы векторизовать этот код, подобно тому, как это сделано в третьем листинге для случая с постоянной

. Но в функции
nchoosek такой возможности нет, приходится все делать в цикле:
n = 6; m = 600;
p = rand(m,n);
p_ = 1-p;
pmf = zeros(m,n+1);
for j=1:m
for i=0:n
P = nchoosek(p(j,:), i);
P_ = flipud( nchoosek(p_(j,:),n-i) );
pmf(j,i+1) = sum( prod([P P_], 2) );
end
end
Для

добавление "небольшой фичи" в код увеличило время его выполнения более чем в два раза. Возможно ли ускорить вычисления? Я так понимаю основная проблема заключается в поиске уникальных комбинаций.
В крайнем случае возможно использование С++.