Discrete probability distribution calculation in Matlab

Question

I have given P(x1...n) discrete independent probability values which represent for example the possibility of happening X.

I want a universal code for the question: With which probability does happening X occur at the same time 0-n times?

For example: Given: 3 probabilities P(A),P(B),P(C) that each car(A,B,C) parks. Question would be: With which probability would no car, one car, two cars, and three cars park?

The answer for example for two cars parking at the same time would be:

P(A,B,~C) = P(A)*P(B)*(1-P(C))
P(A,~B,C) = P(A)*(1-P(B))*P(C)
P(~A,B,C) = (1-P(A))*P(B)*P(C)
P(2 of 3) = P(A,B,~C) + P(A,~B,C) + P(~A,B,C)

I have written the code for all possibilities, but the more values I get, of course the slower it gets due to more possible combinations.

% probability: Vector with probabilities P1, P2, ... PN
% result: Vector with results as stated above.

% All possibilities:
result(1) = prod(probability);


shift_vector = zeros(anzahl_werte,1);
for i = 1:anzahl_werte
    % Shift Vector allocallization
    shift_vector(i) = 1;
    % Compute all unique permutations of the shift_vector
    mult_vectors = uperm(shift_vector);

    % Init Result Vector
    prob_vector = zeros(length(mult_vectors(:,1)), 1);

    % Calc Single Probabilities
    for k = 1:length(mult_vectors(:,1))
        prob_vector(k) = prod(abs(mult_vectors(k,:)'-probability));
    end

    % Sum of this Vector for one probability.
    result(i+1) = sum(prob_vector);
end



end


%%%%% Calculate Permutations
function p = uperm(a)
[u, ~, J] = unique(a);
p = u(up(J, length(a)));
end % uperm

function p = up(J, n)
ktab = histc(J,1:max(J));
l = n;
p = zeros(1, n);
s = 1;
for i=1:length(ktab)
    k = ktab(i);
    c = nchoosek(1:l, k);
    m = size(c,1);
    [t, ~] = find(~p.');
    t = reshape(t, [], s);
    c = t(c,:)';
    s = s*m;
    r = repmat((1:s)',[1 k]);
    q = accumarray([r(:) c(:)], i, [s n]);
    p = repmat(p, [m 1]) + q;
    l = l - k;
end
end
%%%%% Calculate Permutations End

Does anybody know a way to speed up this function? Or maybe Matlab has an implemented function for that?

---

I found the name of the calculation: Poisson binomial distribution

Answer 1

How about this?

probability = [.3 .2 .4 .7];

n = numel(probability);
combs = dec2bin(0:2^n-1).'-'0'; %'// each column is a combination of n values,
%// where each value is either 0 or 1. A 1 value will represent an event
%// that happens; a 0 value will represent an event that doesn't happen.
result = NaN(1,n+1); %// preallocate
for k = 0:n; %// number of events that happen
    ind = sum(combs,1)==k; %// combinations with exactly k 1's
    result(k+1) = sum(prod(...
        bsxfun(@times, probability(:), combs(:,ind)) + ... %// events that happen
        bsxfun(@times, 1-probability(:), ~combs(:,ind)) )); %// don't happen
end