% 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),])