M=80; 
n=10;  % number of balls taken
odds=3; % Prob. of red balls vs other color balls
K=50; % Number of red balls in the urn
samplesize=100000; % Sample size to extract
edges=(0:(n+1))';
% Generate a random number from the distribution above
x=WNChygernd(M,K,n,odds);
disp(x)

M=8000; 
n=2000;  % number of balls taken
odds=3; % Prob. of red balls vs other color balls
K=5000; % Number of red balls in the urn
% Generate an array of size 2x3x5 of random numbers from the distribution above
mm=2;   nn=3; oo=5;
rng(12345)
X=WNChygernd(M,K,n,odds);
disp(X)

M=800; 
n=200;  % number of balls taken
odds=3; % Prob. of red balls vs other color balls
K=500; % Number of red balls in the urn
% Generate an array of size 2x3x5 of random numbers from the distribution above
mm=2;   nn=3; oo=5;
X=WNChygernd(M,K,n,odds, mm,nn,oo);
disp(X)

close all
M=100; % total number of balls
n=10;  % number of balls taken
odds=2; % Prob. of red balls vs other color balls
K=50; % Number of red balls in the urn
samplesize=100000; % Sample size to extract
edges=(0:(n+1))';
% Compute the density
Wpdf=WNChygepdf(edges(1:end-1),M,K,n,odds);
% Generate random numbers
x=WNChygernd(M,K,n,odds,samplesize);
% bar plot of theoretical and relative frequencies
freqANDdens=[(histcounts(x,edges-0.5)/samplesize)' Wpdf];
bar(edges(1:end-1),freqANDdens)
legend(["Theoretical density" "Empirical relative frequency"],'Location','best')
xlabel('Number of successes')