funcalhgpvfishertest.m

calculate the p-value of Fisher Exact Test

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  There a dataset containing N elements, K of which is type A.
%%  randomly select n element form this dataset, i of which is type A.
%%  p-value(k) is the probability that i>=k
%% p-value(k)' is the probability that i=k

% clear

  %%     p12 = P(N21≥k) = sum(Pi)  ,    for i=n : min(N1,N2) ;
  %%      Pi = P(N21=i) = C(N1,i)*C((N-N1),(N2-i))/C(N,N2)
% X=[ 4 7 6 24 ];
% X=[4 44 19 2663 ];
% X=[14 37 142 2562];
% X=[1 3 4 12];
% X=X*1000;
% X=[2 2926 116 62128];

function Pv=funCalHGPvFisherTest(X)
if X(1)>=0 & X(2)>0 & X(2)>=X(1) & X(3)>=X(1) & X(4)>=X(2) & X(4)>=X(3)
    if X(1)==0
        Pv=1;
    else
        k=X(1);n=X(2);K=X(3);N=X(4);
        Pv=0;
        for i=k:min(K,n)
           %%%%%%%  Method 1 (good for its prevention of overflow)
%             A=[1:K,        1:(N-K),            1:n,1:(N-n)];
%             B=[1:i,1:(K-i),1:(n-i),1:(N-K-n+i),1:N];

            A=[(K-i+1):K,   (N-K-n+i+1):(N-K),  (n-i+1):n];
            B=[1:i,         (N-n+1):N      ];
            Pv=Pv+funDivideAs_Bs(A,B);
            
%             Pv=Pv+funHGP(k,n,K,N,i);
           %%%%%%% Method 2
        %    Pv=Pv+nchoosek(K,i)*nchoosek((N-K),(n-i))/nchoosek(N,n);
        end
    % Pv
    end
else
    Pv=NaN;

end

% function Pvp=calHGPvP(i,n,K,N)
% Pvp=nchoosek(K,i)*nchoosek((N-K),(n-i))/nchoosek(N,n);
%