MATLAB: биномиальное распределение с переменной вероятностью

[electric_retard]

Предыстория здесь topic91034.html.
Для понимания вставлю решение с постоянной $p$. Для биномиального распределения с заданными $n$ и $p$ функция вероятности находится следующим образом:

Используется синтаксис Matlab M

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

Теперь просто интегрируем функцию вероятности в требуемом диапазоне. Например, чтобы найти $P(n\geq2)$:

Используется синтаксис Matlab M

n=2;
disp( sum(pmf(n+1:end)) );

Векторизация в этом случае не составляет труда. Если необходимо найти функции распределения сразу для $m$ серий экспериментов:

Используется синтаксис Matlab M

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 производится один раз за цикл, никаких проблем.
В случае переменной вероятности $p$ необходимо для каждого $n_i$ найти все возможные комбинации вероятностей $p$ и $1-p$. На помощь приходит все та же функция nchoosek, а точнее второй вариант ее использования, когда первый аргумент — вектор:

Используется синтаксис Matlab M

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

Первая и вторая строка цикла выдают все возможные комбинации $p$ и $1-p$ в количестве $i$ и $n-i$ раз соответственно. Хотелось бы векторизовать этот код, подобно тому, как это сделано в третьем листинге для случая с постоянной $p$. Но в функции nchoosek такой возможности нет, приходится все делать в цикле:

Используется синтаксис Matlab M

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

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