%  Hypergeometric Probabilities
%  Counting a binary outcome when sampling from
%  finite populations without repetition.
%  Comparison of Theoretical, MATLAB, and Simulation
%  results for P(X=2) when X ~ HG(30, 10, 6)
% needs: histp.m
clear all
close all
  disp('Hypergeometric Example')
  lw = 3; 
  set(0, 'DefaultAxesFontSize', 16);
  fs = 15;
  msize = 10;
 % Notation hygepdf(x, M, K, N) probability of drawing exactly x out of  
 %  K items (x <= N,  x <= K)
 % in N drawings without replacement from a group of M objects.
  M = 30; %  population size
  K = 10; %  number with a particular disease in population
  N = 6;   %  number of people selected at random
  x = 2;   %  number with disease among selected
  
  
 % theoretically
theo_probof2 = nchoosek(K,x) * nchoosek(M-K, N-x)/nchoosek(M, N)
 
 % MATLAB built-in command
matlab_probof2 = hygepdf(2, 30, 10, 6)
  
 % simulation
NN = 10000; %number of simulations
hypers=[];
% construct the population: 1 cases; 0 non-cases
cases = ones(1, K);
nocases = zeros(1, M-K);
popul = [cases nocases];

for i = 1:NN  %repeat the sampling NN times
  sample = randsample(popul, N);
      % if you do not have randsample in matlab, do the following
      % > aa=randperm(N);
      % > selected_indices = aa(1:k);
      % > sample = popul(selected_indices);
      %
  hypers = [hypers sum(sample)];  %find and add x to the vector hypers
end
simul_probof2 = sum(hypers==2)/NN % and estimate the probability of x=2.
%hist(hypers,7)
histp(hypers, 7, 0, 7)
