Poisson binomial distribution MATLAB

Question

could any one please help me to write the "Poisson binomial distribution MATLAB" i have difficulties to do it.

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

Binomial is given as:

function y=mybinomial(n,p)

for k=0:n
    y(k+1)=factorial(n)/(factorial(k)*factorial(n-k))*p^k*(1-p)^(n-k)       
end

i was trying to solve it by this code but it dosnt work Any one could please help me?

clc
clear
k=0;
n=2;
p=[0.7 0.8 ];
prodd=1;
%     pr= zeros(k+1,1);
k=0:n ;
for z=0:length(k);

if k==0;
   for i=1:n  
     prodd=prodd*(1-p(i));
    end

end
pr= prodd

if k > 0 ;

%%%%%%T(i) start  

  for i=1:k;        
    for j=1:n;

        s(j)  =  (p(j)  /   (1-p(j))).^i ;   
    end
    T(i) = sum (s);
  end
%%%%%%%%t(i) end      

 pr= prodd  

 for i=1:length(k)-1; 
  pr(i+1) =   pr(i)+  ((-1)^(i-1))   *     T(i); %%%%%%%%%%%problem        
end    
%     pr(k)=          ((sum(q))/k);

end
end

Many thanks.

Answer 1

Disclaimer: you are trying to implement an approximated formula that only provides acceptable results for a small size of the vector of probabilities.

Your implementation has several problems. You are initializing k as a vector, but it is then used as a scalar. This is also reported by MATLAB as a Warning! This implementation does do what you are asking:

function pmf = poisson_binomial_distribution(p_v)
% approximated, and length(p_v) is small.

inft = length(p_v)+1;
pmf = zeros(1, inft);

% Case K=0
pmf(1) = prod(1-p_v);

% find T(i)
T = zeros(1, inft-1);
for k=1:(inft-1)
    T(k) = sum((p_v ./ ( 1 - p_v )).^(k));
end

% Case k > 0
for k=1:inft-1
    for u=1:k
        pmf(k+1) = pmf(k+1) + (-1)^(u-1)*pmf(k-u+1)*T(u);
    end
    pmf(k+1) = pmf(k+1) / k;
end

pmf(pmf<0)=0;

end

You can see by yourself that this implementation is not always numerically stable, when length(p_v) is large.