Simulated-trials probability not converging to theoretical probability Premise: Person A picks five letters (no repeats) of the standard English alphabet (26 letters). Six other people (persons B-G) are asked to independently pick (privately write on a piece of paper) one letter of the alphabet. Since these picks are independent, it is possible to have repeats in this set.
The task at hand is to calculate the probability that at least one of persons B-G has picked a letter that matches one of the letters picked by person A.
Theoretically, since this probability should be equal to the complement of the probability that none of persons B-G picked a matching letter:
P(at least one match in 6 tries) = 1 - P(no matches in six tries)
                      = 1 - (6C6) * P(not getting a match in a single trial)^6
                      = 1- 1 * (21/26)^6
                      = 0.7224
I set out to generate one million simulated random trials to examine the empirical probability, with the hope that it would come close to the theoretical probability calculated above. My procedure was as follows:


*

*To simulate the action of person A, used MATLAB to generate one million trials of five unique random numbers between 1 to 26, inclusive:


 
for i=1:1000000
    data1(i,:)=randperm(26,5);
end 


*To simulate the actions of persons B-G, I once again used MATLAB to generate one million trials of six random numbers (no unique-number restriction) between 1 to 26, inclusive:



data2 = randi([1 26],1000000,6);



*I quality-checked the data, found that all constraints were satisfied, then exported it to Excel, matching each trial from data set one with a trial from data set two.

*Noting that person A's letter picks are stored in Excel columns A-E, and persons B-G's picks are stored in columns G-L (skipped column F as a divider),  I created a column in Excel that used the following logic to examine each row:

=IF(OR(A1=G1,A1=H1,A1=I1,A1=J1,A1=K1,A1=L1,B1=G1,B1=H1,B1=I1,B1=J1,B1=K1,B1=L1,C1=G1,C1=H1,C1=I1,C1=J1,C1=K1,C1=L1,D1=G1,D1=H1,D1=I1,D1=J1,D1=K1,D1=L1,E1=G1,E1=H1,E1=I1,E1=J1,E1=K1,E1=L1),1,0)

NOTE: above code was edited as per user:lulu's comment below


*I once again quality-checked the results. If one or more of person's B-G's picks matched, the analysis column (column M) displayed a "1". Otherwise, it displayed a zero. Working as intended.

*In order to compute the empirical probability, I performed the following calculation in a separate cell:

=SUM(M:M)/COUNT(M:M)

To my disappointment, the value was 0.67155. I repeated the entire simulation seven more times, and ended up with values of {0.67149, 0.67161, 0.67131, 0.67092, 0.67190, 0.67116, and 0.67126}. This would seem to suggest that the probability is converging to somewhere just above 0.67.
NOTE: with above edit in step 4, for one set of trials, I recomputed to 0.6783. Closer, but still off.
My question is: Is this discrepancy between the theoretical and empirical probabilities due to:


*

*An issue with my simulation methodology in that I incorrectly represented the premises of the initial situation?

*An issue with the simulation methodology in that I did not perform enough trials? (I am doubtful of this one)

*An issue with the pseudorandom number generation?

*An incorrect assessment and computation of the theoretical probability?

*Anything else I have not considered?
I would very much appreciate any suggestions. Thank you!
 A: I'm not sure why you are not doing everything in Matlab (or how you're moving the data over to Excel, which could introduce errors if you're not careful). In pure Matlab:
rng(1); % Set seed to be able to replicate
n = 1e6;
data1 = zeros(n,5);
for i = 1:n
    data1(i,:) = randperm(26,5);
end
data2 = randi([1 26],n,6);

Then, using a double for loop for clarity:
s = 0;
v = zeros(5,1);
for i = 1:n
    for j = 1:5
        v(j) = any(data1(i,j)==data2(i,:));
    end
    s = s+any(v);
end
p = s/n

This returns 0.721587, which is approximately $1-(21/26)^6 = 0.722364062753467...$. Different initial seed will result in different outputs. The double for loop could be vectorized to:
data3 = zeros(n,5);
for i = 1:5
    data3(:,i) = any(bsxfun(@eq,data1(:,i),data2),2);
end
p = sum(any(data3,2))/n

Given this, I'd say the issue is in your Excel code and/or which rows and columns are being used. Did you fill down column M correctly, i.e., do the indices correspond to the row correctly? It's hard to say more without seeing your spreadsheet.
