Расчет дискретного распределения вероятностей в Matlab

Question

Я дал P(x1...n) дискретных независимых значений вероятности, которые представляют, например, возможность возникновения X.

Я хочу универсальный код для вопроса: с какой вероятностью происходит X одновременно в 0-n раз?

Например: Дано: 3 вероятности P(A),P(B),P(C), что каждый автомобиль (A,B,C) паркует. Вопрос будет таким: с какой вероятностью не будет машины, одной машины, двух машин и трех машин?

Например, ответ на парковку двух автомобилей одновременно:

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)

Я написал код для всех возможностей, но чем больше значений я получаю, тем медленнее он становится из-за большего количества возможных комбинаций.

% 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

Кто-нибудь знает способ ускорить эту функцию? Или, может быть, в Matlab есть реализованная функция для этого?

---

Нашел название расчета: биномиальное распределение Пуассона

Answer 1

Как насчет этого?

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