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].