function question()

% Number of yes/no experiments (Bernoulli trials)
N = 79;

% Alpha parameter of the beta distribution
a = 0.85;

% Beta parameter of the beta distribution
b = 0.35;

% Success probability for each of the N trials.
% The p's are not constant, but follow a beta distribution
p = betarnd(a,b,1,N);

% This will store the number of successes
nof_successes = -1*ones(1,200000);

for t = 1:length(nof_successes)
    
    % Generate a uniform random number to decide the outcome of each trial
    X0 = rand(1,N);
    
    % Decide the outcome of each trial (1 = success, 0 = failure)
    S = 1*(X0 <= p);
    
    % Count the number of sucesses in the N trials
    nof_successes(t) = sum(S);
end

figure

% Compute the PDF of k (number of successes) from the simulation
[f,x] = ecdf(nof_successes);
n = ecdfhist(f,x,min(x):max(x));

% Plot the PDF obtained from the simulation and compare to PDF models
plot(min(x):max(x), n, 'b-',...
     0:N, poisson_binomial_PMF(p), 'ro',...
     0:N, beta_binomial_PMF(N, a, b), 'k--')
legend('Simulation', 'Poisson''s binomial', 'Beta-binomial')


%--------------------------------------------------------------------------

function P = poisson_binomial_PMF(p)
%Poisson's binomial probability mass function (PMF)
%P = poisson_binomail_PMF(p)
%   Given N-element vector p with the non-zero probabilities of success of
%   each of N individual, independent trials, this function provides:
%
%   P = N+1 element vector containing the probabilities of 0 successes, 1
%       success, ..., N successess (i.e., the entries of the Poisson-Binomial pdf)

% This implementation has been taken from:
% Fernandez, M.; S. Williams (2010). "Closed-Form Expression for the
% Poisson-Binomial Probability Density Function". IEEE Transactions on
% Aerospace Electronic Systems 46: 803–817. doi:10.1109/TAES.2010.5461658
% (http://ieeexplore.ieee.org/xpl/articleDetails.jsp?arnumber=5461658)

% Eliminate any zero elements of vector p
p(p==0) = [];

N = length(p);

alpha = prod(p);
s = -(1-p)./p;
S = poly(s);            % S = vector with the N+1 cofficients of the polynomial
                        %     whose roots are the elements of vector s; that is,
                        %     S(1)*x^N + ... + S(N)*x + S(N+1)

temp_P = alpha*S;
P = temp_P(N+1:-1:1);   % Reverse the order of the elements of temp_P

%--------------------------------------------------------------------------

% THIS IS THE DISTRIBUTION I'M INTERESTED IN!!!

function P = beta_binomial_PMF(N, alpha_param, beta_param)
% See article in wikipedia
% http://en.wikipedia.org/wiki/Beta-binomial_distribution

% Binomial coefficient is computed as (more efficient, no problems):
% nchoosek(N,k) = round(exp(gammaln(N+1)-gammaln(k+1)-gammaln(N-k+1)))
% http://en.wikipedia.org/wiki/Binomial_coefficient
for k = 0:N
    binomial_coefficient = round(exp(gammaln(N+1)-gammaln(k+1)-gammaln(N-k+1))); % equals nchoosek(N,k)
    P(k+1) = binomial_coefficient*beta(k+alpha_param,N-k+beta_param)/beta(alpha_param,beta_param);
end