function prob=LMTheorem030202_4th(r,k,n,N);
% format prob=LMTheorem030202_4th(r,k,n,N);
% Hypergeometric distribution:
% Suppose an urn contains r red chips and (N-r) white chips
% and n chips are drawn at random and X denotes the total 
% number of red chips selected p_X(k)=P(X=k)=(rCk {N-r}C{n-k}/NCn 
% prob=probability of k successes in n trials
% with probability p.
% input:
%     r=number of red chips in urn
%     N=total number of chips in urn
%     k=scalar or vector, number of red chips withdrawn in sample of size n
%     n=scalar or vector, number of items withdrawn without replacement
%  output:  probability or probability vector
% reference: Theorem 3.3.1, p. 128 in Larsen & Marx (2001)
%            Theorem 3.2.2 p. 139 in Larsen & Marx 4th edition(2006)
% written by E. Gallagher, Environmental Sciences Program
% UMASS/Boston, written: 5/1/03, last revised: 10/3/2010
% Email: Eugene.Gallagher@umb.edu
% Note that the Matlab statistics toolbox has hygepdf.m, which uses
% Stirling's approximation to the factorial instead of gammaln

if k>r
    error ('r < k')
elseif n>N
    error('N < n')
end
flipt=0;
rows=1;
if length(n)>1 & length(k)>1
    error('k or n must be a scalar')
elseif length(k)>1
    [rows,cols]=size(k);
    if rows<cols
        flipt=1;
        k=k(:);
        rows=cols;
    end
    n=repmat(n,rows,1);
elseif length(n)>1
    [rows,cols]=size(n);
    if rows<cols
        flipt=1;
        n=n(:);
        rows=cols;
    end
    k=repmat(k,rows,1);
end
w=N-r;  % White chips = total (N) - red;
w=repmat(w,rows,1);
N=repmat(N,rows,1);
r=repmat(r,rows,1);
onek=ones(rows,1);
prob=zeros(rows,1);
i=find((w-(n-k)+onek)>0);j=find((w-(n-k)+onek)<=0);
if length(j)>0;prob(j)=0;end
onek=ones(length(i),1);
lnNf=gammaln(N(i)+onek);
lnwf=gammaln(w(i)+onek);
lnwmnmkf=gammaln(w(i)-(n(i)-k(i))+onek);
lnrf=gammaln(r(i)+onek);
lnkf=gammaln(k(i)+onek);
lnnf=gammaln(n(i)+onek);
lnNmnf=gammaln(N(i)-n(i)+onek);
lnrmkf=gammaln(r(i)-k(i)+onek);
lnnmkf=gammaln(n(i)-k(i)+onek);
prob(i)=exp( (lnrf-(lnkf+lnrmkf) + lnwf-(lnnmkf+lnwmnmkf) )- ...
    (lnNf -(lnnf+lnNmnf)));
if flipt==1;prob=prob';end