nhygepdf(x,m,n,a)

function y = nhygepdf(x,m,n,a)
%NHYGEPDF Negative hypergeometric probability density function.
% Y = NHYGEPDF(X,M,N,A) returns the negative hypergeometric probability
% density function with parameters M, N and A at the values in X. 
% Note: The density function is zero unless M, N and A are integers.
%
% If a lot consists of M acceptable items and N defective ones. Suppose
% that items are drawn at random one by one without replacement and x 
% defectives are observed before the a'th acceptable one. According to 
% Johnson et al. (2005), in an inspection sampling, instead of taking a
% sample of fixed size from a batch of items and then accepting the batch
% if the observed number of defectives is less than or equal to some 
% predetermined value (otherwise rejecting). Here, items are drawn one at a
% time until either defectives are observed (at which point the batch is
% rejected) or nondefectives are observed (and the batch is accepted). The 
% number of observations required to reach a decision then has a negative
% hypergeometric distribution.
%
% This discrete distribution is also known as Beta-Binomial, Inverse
% hypergeometric, Hypergeometric waiting time, and Markov-Polya.
%
% Its distribution function is expresed as,
%
%          f(x) = P(X=x) = (x+a-1_C_x)(m+n-a-x_C_m-x)/(m+n_C_m)
%
%     'Probability of x failures one will have before the a success(es)'
%
% where: m=number of acceptables, n=number of not acceptables (defectives),
% x=number of observed defectives, before a=acceptable ones. C means 
% combination. m+n=N, total number of elements.
%
% Syntax: function y = nhygepdf(x,m,n,a) 
%      
% Inputs:
% m - number of acceptables  
% n - number of not acceptables (defectives)  
% x - number of observed defectives, before
% a - acceptable ones
%
% Output:
% y - negative hypergeometric probability value
%
% Example 1. From the example given on the web page PlanetMath 
%  [http://planetmath.org/encyclopedia/ExampleOfNegativeHypergeometricRandomVariable.html].
%  Suppose have 7 black marbles and 10 white marbles in a jar. You pull 
%  marbles until you have 3 black marbles in your hand. X would represent
%  the number of white marbles in your hand. So, we are interested to 
%  estimate the probability of having 3 white marbles before of 3 black
%  marbles.   
%
%              x=3; m=10; n=7; a=3;
%
% Calling on Matlab the function: 
%             y = nhygepdf(x,m,n,a)
%
% Answer is: (in format long)
%
% y =
%
%   0.16968325791855
%
% Created by A. Trujillo-Ortiz, R. Hernandez-Walls, F.A. Trujillo-Perez
%            and N. Castro-Castro
%            Facultad de Ciencias Marinas
%            Universidad Autonoma de Baja California
%            Apdo. Postal 453
%            Ensenada, Baja California
%            Mexico.
%            [e-mail]
% Copyright. September 25, 2009.
%
% ---Con cario para Ney.----
%
% To cite this file, this would be an appropriate format:
% Trujillo-Ortiz, A., R. Hernandez-Walls, F.A. Trujillo-Perez and 
%   N. Castro-Castro (2009). nhygepdf:Negative hypergeometric 
%   probability density function. A MATLAB file. [WWW document].
%   URL http://www.mathworks.com/matlabcentral/fileexchange/25434
%
% Reference:
% Johnson, N.L., A.W. Kemp, and S. Kotz. (2005), Univariate Discrete
%      Distributions. Wiley Interscience. John Wiley & Sons:New York.
%


if nargin < 4, 
    error('nhygepdf:TooFewInputs','Requires four input arguments.'); 
end

[errorcode x m n a] = distchck(4,x,m,n,a);

if errorcode > 0
    error('nhygepdf:InputSizeMismatch',...
        'Requires non-scalar arguments to match in size.');
end

if (length(x)~=1) || (fix(x) ~= x) || (x < 0),
   error('nhygepdf:InvalidData',...
       'NHYGEPDF requires that X must be a non-negative and integer.')
end

if (length(m)~=1) || (fix(m) ~= m) || (m < 0),
   error('nhygepdf:InvalidData',...
       'NHYGEPDF requires that M must be a non-negative and integer.')
end

if (length(n)~=1) || (fix(n) ~= n) || (n < 0),
   error('nhygepdf:InvalidData',...
       'NHYGEPDF requires that N must be a non-negative and integer.')
end

if (length(a)~=1) || (fix(a) ~= a) || (a < 0),
   error('nhygepdf:InvalidData',...
       'NHYGEPDF requires that A must be a non-negative and integer.')
end

%y = exp(gammaln(x + a -1 + 1) - gammaln(x + 1) - gammaln(x + a -1 - x + 1) +...
%    gammaln(m + n - a - x + 1) - gammaln(m - x + 1) - gammaln(m + n - a - x - m + x + 1) -...
%    (gammaln(m + n + 1) - gammaln(m + 1) - gammaln(m + n - m + 1)));

%After an algebraic simplification we have
N = m + n;

y = exp(gammaln(x + a) + gammaln(N - a - x + 1) + gammaln(m + 1) +...
    gammaln(n + 1) - gammaln(x + 1) - gammaln(a) - gammaln(m - x + 1) -...
    gammaln(n - a + 1) - gammaln(N + 1));

return,