MATLAB Answers

Simulation of beta-binomial distribution

28 views (last 30 days)
Hi all!
I'm trying to solve the following problem.
The number of successes in a sequence of N yes/no experiments (i.e., N Bernoulli trials), each of which yields success with probability p, is given by the well-known binomial distribution. This is true if the success probability p is constant and the same for all the N trials.
However, when the probability of success, p, is different for each trial, p_1, p_2, ..., p_N, then the number of successes does not follow a binomial distribution, but a Poisson's binomial distribution instead:
I understand that the Poisson's binomial distribution is valid for any set of probabilities p_1, p_2, ..., p_N.
In the problem I'm trying to solve, the probabilities p_1, p_2, ..., p_N follow a beta distribution and I know their values. Can this knowledge be exploited to obtain a PMF that is simpler than the Poisson's binomial PMF?
I found out that the PMF of the number of successes in N trials where the success probability is a beta-distributed random variable is given by the beta-binomial distribution:
However, I have been playing a bit with some simulation and it seems that this distribution does not fit the resulting PMF. I'm attaching below the Matlab code that makes some simulation and generates the PMFs.
What am I doing wrong? Is it possible to exploit the knowledge that the p_1, p_2, ..., p_N follow a beta distribution to simplify the general Poison's binomial case? What is the PMF that I need?
Many thanks in advance!
Fryderyk C.
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);
% 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 % (<>)
% 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 coefficients 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
function P = beta_binomial_PMF(N, alpha_param, beta_param) % See article in wikipedia %
% Binomial coefficient is computed as (more efficient, no problems): % nchoosek(N,k) = round(exp(gammaln(N+1)-gammaln(k+1)-gammaln(N-k+1))) %
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);

Accepted Answer

Lokesh Bommisetty
Lokesh Bommisetty on 21 Oct 2019
The following implementation will give the pmf of poisson binomial distribution accuratly with the lest complexiety.
function S=Poisson_Binomial(p,N)
for k=0:N
for l=0:N

More Answers (1)

Paul Fackler
Paul Fackler on 27 Jan 2014
I believe that the distribution you are looking for is simply Binomial(n,a/(a+b))

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by