% LMex030207_4th.m
% Example 3.2.7, Bullet striations on page 142
% Larsen & Marx (2006) Introduction to Mathematical Statistics, 4th edition
% Written by Eugene.Gallagher@umb.edu, written in 2001, revised 1/12/11
% http://alpha.es.umb.edu/faculty/edg/files/edgwebp.html
% Solution to bullet striation pattern problem
% will call Matlab statistics toolbox hypergeometric probability
% distribution function or Gallagher's hypergeop.m
% The problem will also be solved using a Monte Carlo simulation.
% There were 4 striations in 25 locations in the evidence bullet
% A test bullet from the suspect's gun was examined and 1 striation out of
% 3 on the test bullet matched the striations on the evidence bullet.
% What is the probability that 1 or or more of the test bullet striations
% would match the test bullet
if exist('hygepdf','file')
    P=sum(hygepdf(1:3,25,3,4));
elseif exist('hypergeop','file')
    P=sum(hypergeop(25,3,4,1:3));
else
    disp('No hypergeometric function on path')
end
fprintf(...
'Probability of 1, 2 or 3 matches from hypergeometric distribution =%7.5f.\n',...
   P)
% According to the book, small p values indicate that the evidence
% gun and the test gun were the same.
% What if 2 of the 3 patterns from the test bullet matched the
% evidence bullet?
if exist('hygepdf','file')
   NewP=sum(hygepdf(2:3,25,3,4));
elseif exist('hypergeop','file')
   NewP=sum(hypergeop(25,3,4,2:3));
else
    disp('No hypergeometric function on path')
end
fprintf(...
'Probability of  2 or 3 matches from hypergeometric distribution =%7.5f.\n',...
   NewP)
% What if all 3 of the 3 patterns from the test bullet matched the
% evidence bullet?
if exist('hygepdf','file')
   NewP=sum(hygepdf(3,25,3,4));
elseif exist('hypergeop','file')
   NewP=sum(hypergeop(25,3,4,3));
else
    disp('No hypergeometric function on path')
end
fprintf(...
'Probability of 3 matches from hypergeometric distribution =%7.5f.\n',...
   NewP)
% Advanced problem:
% This problem is amenable to a Monte Carlo simulation
EB=[ones(1,4) zeros(1,21)]; % evidence bullet with 4 striations
TB=[1 0 0 0 1 1 zeros(1,19)];% test bullet with 3 striations, 1 overlapping
trials=9999;
matches=sum((EB+TB)==2);
teststat=sum(hygepdf(matches:3,25,3,4));
results=zeros(trials,1);
for i=1:trials
    % just need to shuffle one of the bullets, just the evidence bullet
    EBshuffled=EB(randperm(length(EB)));
    %    results(i)=sum((EBshuffled+TB)==2); % Probability of 1 match
    results(i)=sum((EBshuffled+TB)>=2);  % What is the probability of 1 OR
                                         % MORE matches
end
i=find(results>=matches);MCP=(length(i)+1)/trials;
fprintf(...
  'P from hypergeometric =%7.5f & from Monte Carlo = %6.4f +/- %6.4f\n',...
   P, MCP,norminv(0.975)*sqrt(P*(1-P))/sqrt(trials))
% Note that if there HAD been more than 1 match, the spacing of the 
% striations would play a key role. The random permutations would have
% to be performed taking into account the spatial autocorrelation of
% striations on the bullets. ter Braak discusses how to shuffle in the
% presence of spatial autocorrelation.