Non-Uniform Bernoulli Trials Exceedance Probability

Calculates the probability of exceeding a given number of hits for a collection of Bernoulli trials

Description

Calculates the probability of exceeding a given number of hits for a collection of Bernoulli trials, which is a test with a binary result (e.g. yes/no, 1/0, heads/tails). User inputs a vector of probabilities of success for each trial, p, and a number of successes, h, and this function outputs the probability of equaling or exceeding h.

Example:
There is a 30% chance of rain Saturday, and 20% Sunday. What is the
probability that it rains this weekend?

monte_exceedance([.3; .2],1)

ans =
    0.4431

Which is reasonably close to the theoretical value of .44. Increasing the number of runs improves the accuracy:

monte_exceedance([.3; .2],1,100000)

ans =
    0.4407

• MATLAB release: MATLAB 7.5 (R2007b)

[Darren Rowland]

The entire contents of the file...
Frankly the quality is not worthy of the FEX. Please remove.

function exc = monte_exceedance( p, h, runs )
%
% MONTE_EXCEEDANCE( P, H, RUNS ) calculates the probability of equaling or
% exceeding h hits for probabilities p (column vector). Uses a monte carlo method,
% as the author finds analytical statistics tiresome. This should reflect
% binomial distribution calculations for uniform p.
%
% Example:
% There is a 30% chance of rain Saturday, and 20% Sunday. What is the
% probability that it rains this weekend?
%
% monte_exceedance([.3; .2],1)
%
% ans =
% 0.4431
%
% Which is reasonably close to the theoretical value of .44

if nargin < 3
    runs = 10000;
end

r = rand(length(p),runs);
exc = mean( sum(bsxfun(@lt,r,p)) >= h );

end

[William]

That is rather frank...
Could you point me to a better way to perform this calculation, perhaps an analytical method?
Thanks.

[Peter Cotton]

William,

   I was looking for the same thing and ended up re-writing it myself (bernpdf.m). Let me know if you come across other analytic implementations (e.g. by Fourier transform).

   I think you mean 'non-uniform p', for otherwise you can use binopdf.m in the stats toolbox.

Peter

ps: Sorry you find analytical statistics tiresome