pdf of Poisson binomial distribution in Matlab

Is there a Matlab implementation for the calculation of the pdf for the Poisson-binomial distribution?

[Torsten]

No. How large is n in your case ?

The usual recusion to compute

P( K = k ) for 0 <= k <= n

become unstable for n large:

https://en.wikipedia.org/wiki/Poisson_binomial_distribution

[valentino dardanoni]

Thank you Torsten, I guess I could start with a small example with small n, say n=5

[Torsten]

There are some MATHEMATICA codes that you could compare your results against:

https://mathematica.stackexchange.com/questions/196962/does-mathematica-have-an-implementation-of-the-poisson-binomial-distribution

Answer 1

[Gyan Vaibhav]

Hi Valentino,

I understand that you are trying to find the PDF of a given Binomial-Poisson distribution.

While MATLAB doesn't offer a built-in function specifically for this purpose, you can certainly craft a custom function to accomplish the task.

The code snippet provided below is designed to calculate the PDF for a Poisson-Binomial distribution. This function requires two input arguments:

1. successProbs: A vector containing the individual success probabilities for each trial.
2. k: The specific number of successful trials for which you wish to compute the PDF.

function pdf = poisson_binomial_pdf(successProbs, k)
% successProbs is a vector containing the success probabilities for each trial
% k is the number of successful trials for which you want to calculate the PDF

n = length(successProbs); % Number of trials
% The FFT-based method for Poisson-binomial PDF calculation
M = 2^nextpow2(2*n); % Find the next power of 2 for zero-padding
omega = exp(-2i * pi / M);
A = ones(1, M);
for j = 1:n
    A = A .* (1 - successProbs(j) + successProbs(j) * omega.^(0:M-1));
end
pdfVals = ifft(A);
pdf = real(pdfVals(1:n+1)); % Only the first (n+1) values are needed

% Return the PDF value for k successes
if k >= 0 && k <= n
    pdf = pdf(k+1);
else
    pdf = 0;
end
end

You can use this function as follows:

% Example success probabilities for 5 trials
successProbs = [0.04, 0.07, 0.07];
% Calculate the PDF for 3 successes
k = 3;
pdfValue = poisson_binomial_pdf(successProbs, k);

This approach gives us the PDF of a Binomial-Poisson Distribution.

Hope this helps.

Thanks

Gyan

[Paul]

Hi Gyan,

Thanks for posting your function. I ran a quick test to see it work.

Question: If I understand correctly, could the function just return after this line

pdf = real(pdfVals(1:n+1));

and pdf would be the probabilities for all values of k in 0 <= k <= n ?

Also, it may be useful to consider handling the special case where one or more elements of successProbs are exactly zero.

rng(100);
successProbs = rand(1,10);
ntrials = numel(successProbs);
for k = 0:ntrials
    pdf(k+1) = poisson_binomial_pdf(successProbs,k);
end
N = 1e5;
trials = rand(N,ntrials);
success = trials < successProbs;
for k = 0:ntrials
    epdf(k+1) = sum(sum(success,2) == k)/N;
end
figure
plot(0:ntrials,pdf,'-x',0:ntrials,epdf,'-o')

function pdf = poisson_binomial_pdf(successProbs, k)
% successProbs is a vector containing the success probabilities for each trial
% k is the number of successful trials for which you want to calculate the PDF
n = length(successProbs); % Number of trials
% The FFT-based method for Poisson-binomial PDF calculation
M = 2^nextpow2(2*n); % Find the next power of 2 for zero-padding
omega = exp(-2i * pi / M);
A = ones(1, M);
for j = 1:n
    A = A .* (1 - successProbs(j) + successProbs(j) * omega.^(0:M-1));
end
pdfVals = ifft(A);
pdf = real(pdfVals(1:n+1)); % Only the first (n+1) values are needed
% Return the PDF value for k successes
if k >= 0 && k <= n
    pdf = pdf(k+1);
else
    pdf = 0;
end
end

[valentino dardanoni]

Thank you Gyan and Paul

[Gyan Vaibhav]

Hi Paul,

Thanks for the cross-check.

Regarding your question, yes you are correct about that. If we remove the last if-else statement we can get the probabilities for k in the range [0,n].

[John D'Errico]

Note, this code was at least viable for n on the order of 1e4. But for 1e5 (the case I wanted to use it for, thus with vectors of binomial probabilites that were of length 1e5) it is pretty slow.

[Paul]

Hi John,

I was able to get significant speedup in the code by hoisting the omega.^ out of the loop, using the 'symmetric' flag, and using the fact that we don't need to loop over k. In addition, I put comment in the new function that shows what's probably a better way to compute the omega.^ term, but left it commented so as to make sure I got the exact same results as the original function.

I did not try the code for the scope of problem that you're interested in.

I wonder if further gains can be made in the computation of A.

rng(100);
successProbs = rand(1,500);
ntrials = numel(successProbs);
pdf1 = nan(ntrials+1,1);

Original code

tic
for k = 0:ntrials
    pdf1(k+1) = poisson_binomial_pdf(successProbs,k);
end
toc

Elapsed time is 10.251024 seconds.

New code

tic
pdf2 = poisson_binomial_pdfnew(successProbs,ntrials);
toc

Elapsed time is 0.008475 seconds.

Estimated result

tic
N = 1e5;
trials = rand(N,ntrials);
success = trials < successProbs;
for k = 0:ntrials
    epdf(k+1) = sum(sum(success,2) == k)/N;
end
toc

Elapsed time is 2.646702 seconds.

Compare results

max(abs(pdf1(:)-pdf2(:)))

ans = 0

figure
plot(0:ntrials,pdf1,'-x',0:ntrials,pdf2,'-o',0:ntrials,epdf,'-d')
legend('pdf1','pdf2','epdf')

% original function
function pdf = poisson_binomial_pdf(successProbs, k)
% successProbs is a vector containing the success probabilities for each trial
% k is the number of successful trials for which you want to calculate the PDF
n = length(successProbs); % Number of trials
% The FFT-based method for Poisson-binomial PDF calculation
M = 2^nextpow2(2*n); % Find the next power of 2 for zero-padding
omega = exp(-2i * pi / M);
A = ones(1, M);
for j = 1:n
    A = A .* (1 - successProbs(j) + successProbs(j) * omega.^(0:M-1));
end
pdfVals = ifft(A,'symmetric');
pdf = (pdfVals(1:n+1)); % Only the first (n+1) values are needed
% Return the PDF value for k successes
if k >= 0 && k <= n
    pdf = pdf(k+1);
else
    pdf = 0;
end
end

% new function
function pdf = poisson_binomial_pdfnew(successProbs, k)
% successProbs is a vector containing the success probabilities for each trial
% k is the number of successful trials for which you want to calculate the PDF
n = length(successProbs); % Number of trials
% The FFT-based method for Poisson-binomial PDF calculation

% should this be nextpow2(n)?
M = 2^nextpow2(2*n); % Find the next power of 2 for zero-padding
omega = exp(-2i * pi / M);
omegaM = omega.^(0:M-1);
% maybe would be better to do
% omegaM = exp(-2i*pi*(0:M-1)/M);
A = ones(1, M);
for j = 1:n
    A = A .* (1 - successProbs(j) + successProbs(j) * omegaM);
end
pdfVals = ifft(A,'symmetric');
pdf = pdfVals(1:n+1); % Only the first (n+1) values are needed

%{
% Return the PDF value for k successes
if k >= 0 && k <= n
    pdf = pdf(k+1);
else
    pdf = 0;
end
%}

end

[John D'Errico]

@Paul

Yeah, I did not look carefully at the original code. Pulling omega out helps a lot. Given that I have one case where n may be as large as 7 million, the difference would be significant.

[Paul]

I didn't test each change individually, but I imagine elminating the loop over k (outside the function) is quite signficant as well.