MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn moreOpportunities for recent engineering grads.

Apply Today**New to MATLAB?**

Asked by Fryderyk
on 4 May 2012

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:

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

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:

http://en.wikipedia.org/wiki/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);

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 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

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

% 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

*No products are associated with this question.*

Answer by Paul Fackler
on 27 Jan 2014

I believe that the distribution you are looking for is simply Binomial(n,a/(a+b))

see: http://repositories.lib.utexas.edu/handle/2152/ETD-UT-2012-05-5555

## 0 Comments