binomial distribution with variable probability: code optimisation

Hi

Theory is here. There is no problem to get pmf:

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

Now you have pmf, just sum it over desierable interval to obtain probability that you need. For example: to get probability when n=>2:

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

This case could be also vectorized for multiple p since you need only one nchoosek call for each iteration:

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

I have code which calculates pmf for binomial distribution, where p varies from experiment to experiment:

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

As you can see I need to generate all combinations of p's for each n_i. My aim is to vectorize this code for case when you have m x n different probabilities. Here is my solution for now:

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

So the question is if there is more efficient and gentle way to perform computations from last listing? As i understand I can not vectorize nchoosek call. So it looks like I need vectorized "find all combinations" method. Any ideas?