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.
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 );
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
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
Comment only