pdf of Poisson binomial distribution in Matlab

Question

valentino dardanoni 2023-12-9

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2058364-pdf-of-poisson-binomial-distribution-in-matlab

评论： Paul 2024-12-22

采纳的回答： Gyan Vaibhav

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

3 个评论
显示 1更早的评论隐藏 1更早的评论

valentino dardanoni 2023-12-9

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

Torsten 2023-12-9

编辑：Torsten 2023-12-9

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 2023-12-18

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2058364-pdf-of-poisson-binomial-distribution-in-matlab#answer_1373412

编辑：Gyan Vaibhav 2023-12-18

在 MATLAB Online 中打开

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:

successProbs: A vector containing the individual success probabilities for each trial.
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

7 个评论
显示 5更早的评论隐藏 5更早的评论

Paul 2023-12-18

在 MATLAB Online 中打开

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

Paul 2024-12-22

在 MATLAB Online 中打开

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

% 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 2024-12-22

编辑：John D'Errico 2024-12-22

@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 2024-12-22

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

请先登录，再进行评论。

pdf of Poisson binomial distribution in Matlab

3 个评论
显示 1更早的评论隐藏 1更早的评论

采纳的回答

7 个评论
显示 5更早的评论隐藏 5更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

pdf of Poisson binomial distribution in Matlab

3 个评论 显示 1更早的评论隐藏 1更早的评论

采纳的回答

7 个评论 显示 5更早的评论隐藏 5更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

3 个评论
显示 1更早的评论隐藏 1更早的评论

7 个评论
显示 5更早的评论隐藏 5更早的评论