% CARDMATCH file cardmatch.m Prob of matches in cards from identical decks
% Version of 6/27/97
% Estimates the probability of one or more matches
% in drawing cards from nd decks of c cards each
% Produces a supersample of size n = nd*ns, where
% ns is the number of samples
% Each sample is sorted, and then tested for differences
% between adjacent elements.  Matches are indicated by
% zero differences between adjacent elements in sorted sample
c  = input('Enter the number c  of cards in a deck ');
nd = input('Enter the number nd of decks ');
ns = input('Enter the number ns of sample runs ');
X  = 1:c;                   % Population values
PX = (1/c)*ones(1,c);       % Population probabilities 
N  = nd*ns;                 % Length of supersample
U  = rand(1,N);             % Matrix of n random numbers
T  = dquant(X,PX,U);        % Supersample obtained with quantile function;
                            % the function dquant determines quantile
                            % function values of random number sequence U
ex = sum(T)/N;              % Sample average
EX = dot(X,PX);             % Population mean
vx = sum(T.^2)/N - ex^2;    % Sample variance
VX = dot(X.^2,PX) - EX^2;   % Population variance
A  = reshape(T,nd,ns);      % Chops supersample into ns samples of size nd
DS = diff(sort(A));         % Sorts each sample
m  = sum(DS==0)>0;          % Differences between elements in each sample
                            % Zero difference iff there is a match
pm = sum(m)/ns;             % Fraction of samples with one or more matches
Pm = 1 - comb(c,nd)*gamma(nd + 1)/c^(nd);  % Theoretical probability of match
disp('The sample is in column vector T')   % Displays of results
disp(['Sample average ex = ', num2str(ex),])
disp(['Population mean E(X) = ',num2str(EX),])
disp(['Sample variance vx = ',num2str(vx),])
disp(['Population variance V(X) = ',num2str(VX),])
disp(['Fraction of samples with one or more matches   pm = ', num2str(pm),])
disp(['Probability of one or more matches in a sample Pm = ', num2str(Pm),])