% TRIALMATCH file trialmatch.m  Estimates probability of matches 
% in n independent trials from identical distributions
% Version of 8/20/97
% Estimates the probability of one or more matches 
% in a random selection from n identical distributions
% with a small number of possible values
% Produces a supersample of size N = n*ns, where
% ns is the number of samples.  Samples are separated.
% Each sample is sorted, and then tested for differences
% between adjacent elements.  Matches are indicated by
% zero differences between adjacent elements in sorted sample.
X  = input('Enter the VALUES in the distribution ');
PX = input('Enter the PROBABILITIES  ');
c  = length(X);
n  = input('Enter the SAMPLE SIZE n ');
ns = input('Enter the number ns of sample runs ');
N  = n*ns;                  % Length of supersample
U  = rand(1,N);             % Vector of N random numbers
T  = dquant(X,PX,U);        % Supersample obtained with quantile function;
                            %   the function dquant determines quantile
                            %   function values for 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,n,ns);       % Chops supersample into ns samples of size n
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
d  = arrep(c,n);
p  = PX(d);
p  = reshape(p,size(d));    % This step not needed in version 5.1
ds = diff(sort(d))==0;
mm = sum(ds)>0;
m0 = find(1-mm);
pm0 = p(:,m0);              % Probabilities for arrangements with no matches
P0 = sum(prod(pm0));      
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(1-P0),])