| << Chapter < Page | Chapter >> Page > |
trialmatch.m Estimates the probability of matches in n independent trials from identical distributions. The sample size and number of trials mustbe kept relateively small to avoid exceeding available memory.
% 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 numbersT = 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 averageEX = dot(X,PX); % Population mean
vx = sum(T.^2)/N - ex^2; % Sample varianceVX = dot(X.^2,PX) - EX^2; % Population variance
A = reshape(T,n,ns); % Chops supersample into ns samples of size nDS = diff(sort(A)); % Sorts each sample
m = sum(DS==0)>0; % Differences between elements in each sample
% -- Zero difference iff there is a matchpm = 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.1ds = 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),])
comb.m
function y = comb(n,k) Calculates binomial coefficients.
k may be a matrix of integers between 0 and
n . The result
y is a matrix of the same
dimensions.
function y = comb(n,k)
% COMB y = comb(n,k) Binomial coefficients% Version of 12/10/92
% Computes binomial coefficients C(n,k)% k may be a matrix of integers between 0 and n
% result y is a matrix of the same dimensionsy = round(gamma(n+1)./(gamma(k + 1).*gamma(n + 1 - k))); ibinom.m Binomial distribution — individual terms. We have two m-functions ibinom and cbinom for calculating individual and cumulative terms, and , respectively.
function y = ibinom(n,p,k) calculates the probabilities
for
n a positive integer,
k a matrix of integers between 0 and
n . The output is a matrix of the corresponding binomial probabilities. function y = ibinom(n,p,k)
% IBINOM y = ibinom(n,p,k) Individual binomial probabilities% Version of 10/5/93
% n is a positive integer; p is a probability% k a matrix of integers between 0 and n
% y = P(X>=k) (a matrix of probabilities)
if p>0.5
a = [1 ((1-p)/p)*ones(1,n)];
b = [1 n:-1:1];
c = [1 1:n];
br = (p^n)*cumprod(a.*b./c);bi = fliplr(br);
elsea = [1 (p/(1-p))*ones(1,n)];b = [1 n:-1:1];c = [1 1:n];bi = ((1-p)^n)*cumprod(a.*b./c);
endy = bi(k+1); 
Notification Switch
Would you like to follow the 'Applied probability' conversation and receive update notifications?
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|