# What is an approximation for Poisson binomial distribution?

I am looking for an approximation for Poisson binomial distribution:

The Poisson binomial distribution is the discrete probability distribution of a sum of $n$ independent Bernoulli trials. you can find its pdf in http://en.wikipedia.org/wiki/Poisson_binomial_distribution

In addition, you can find 2 methods for it in the mentioned link. But when I use the second method (using Fourier Transform), the result would be an imaginary number.

I also wanted to use approximations in the following paper http://statistics.stanford.edu/~ckirby/techreports/ONR/SOL%20ONR%20467.pdf

but the approximations are not clear to me. I would be grateful if somebody explains to me:

1) why do I get imaginary number using the second method ( Fourier transform)?
2) and also for example, how are the probabilities in Table 2 in the paper calculated?

-
 "In addition, you can find 2 approximations for it in the mentioned link." I don't see approximations in the wikipedia link, those are exact formulas. – leonbloy Dec 11 '12 at 1:55 @leonbloy: Yes you are right! I edited the question. – May Dec 11 '12 at 15:44 How are we supposed to tell you why you get an imaginary number if you don't tell us how you get it?! – joriki Dec 11 '12 at 16:05 @joriki: For example consider these success probabilities: p1=0.5, p2=0.3, p3=0.7, p4=0.3, p5=0.9. And now we want to calculate P(K=2) . Using The first and second(recursive one) formula (in en.wikipedia.org/wiki/Poisson_binomial_distribution) I get 0.3167, but using the third one (Fourier Transform) I get 0.1957 + 0.0944i. – May Dec 11 '12 at 16:32 As I said, I don't understand how you think we can tell you why you got a wrong result if you don't show us how you got it. The Wikipedia formula is manifestly real-valued, since the contributions for $l$ and $n+1-l$ are complex conjugates of each other, and the remaining contributions for $l=0$, and if $n$ is odd for $l=(n+1)/2$, are real-valued. – joriki Dec 11 '12 at 17:48

1) Post the script/code/etc. you used, there is an error in it that perhaps readers can help you with: I whipped up the example in my CAS and it gets proper answer within the precision set.

2) The examples in that paper use multiples for each Pn, i.e., in its first example with comparisons in table 2, we have Pn of [.02,.04,.06,.08,.10], but there are 5 each accounted for, so you need to apply the calculation to the list of probabilities consisting of five repetitions of each value. Same with later examples. The result tables then show the true value and the various approximations that S(actual expected successes)<= s for some number of s successes, so you need to sum from 0 to s in your calculations to compare your results.

As previous posters stated, the Wikipedia DFT form is not an approximation, it is exact (modulo the capability of your computing system). It will take some serious cycles with large arguments - the form lends itself to using DFT tricks to speed calculations, but if you just implement it "as is" in say a CAS, you're not gaining much efficiency if any over the recursive form sans call stack space.

-
 I posted my code thank you. – May Feb 13 at 16:44

This is My MATLAB code for Fourier Transform method (algorithm A in the paper "On computing the distribution function for the Poisson binomial distribution") and I do not know where the problem is:

function [result,p]=Poisson_binomial_Probability_density_function_FourierTeansform(k,success_probailities)
format long
n=length(success_probailities);
if any(success_probailities<0)|| any(success_probailities>1)
error('invalid values in success_probailities array!');
end
if k>length(success_probailities)
error('invalid value for k!! (0<=k<=length(success_probailities))');
end
a=zeros(1,m);
b=zeros(1,m);
for g=1:m
[a(g),b(g)]=x_real_imaginary(success_probailities,g);
a(n+1-g)=a(g);
b(n+1-g)=-b(g);
x(g)=a(g)+1i*b(g);
x(n+1-g)=a(n+1-g)+1i*b(n+1-g);
end
x=[1,x];
p=fft(x/(n+1));
result=p(k+1);
end

function [result]=abs_z(p,l,n)
%p is the success probability
w=(2*pi)/(n+1);
r1=(1-p+(p*cos(l*w)))^2;
r2=(p*sin(l*w))^2;
result=sqrt(r1+r2);
end

function [result]=Arg_z(p,l,n)
%p is the success probability
r1=1-p+(p*cos(l*(2*pi)/(n+1)));
r2=p*sin(l*(2*pi)/(n+1));
result=atan22(r2,r1);
end

function [result]=atan22(y,x)
if x>0
result=atan(y/x);
elseif y>=0 && x<0
result=atan(y/x)+pi;
elseif y<0 && x<0
result=atan(y/x)-pi;
elseif y>0 && x==0
result=pi/2;
elseif y<0 && x==0
result=-pi/2;
elseif y==0 && x==0
result=0;
end
end

function [result]=d(success_probailities,l)
sum=0;
n=length(success_probailities);
for i=1:n
sum=sum+log(abs_z(success_probailities(i),l,n));
end
result=exp(sum);
end

function [a,b]=x_real_imaginary(success_probailities,l)
r=0;
n=length(success_probailities);
for i=1:n
r=r+Arg_z(success_probailities(i),l,n);
end
dl=d(success_probailities,l);
a=dl*cos(r);
b=dl*sin(r);
end

-