Monte Carlo Method Fails

I am trying to find the probability of getting a pair when drawing 2 cards from a regular deck that contains 52 cards by using Monte-Carlo Method. The function needs one input, a number of trials. However, I noticed that as the input increases the probability converges to zero. I think the theoretical probability of getting a pair is 3/51. Any ideas why my function fails?

function onepair(k)


p = zeros(1,k);

for i = 1:k
    %% Making a regular deck with four patterns and 13 numbers (including A,J,Q,K)
    % randperm is a function that creats an array with N elements that are
    % randomly selected from 1:N
   
    hearts = randperm(13);
    diamonds = randperm(13);
    spades = randperm(13);
    clovers = randperm(13);
    allcards = [hearts,diamonds,spades,clovers];
    
    
    %% Picking up a card
    % datasample is a function that picks random sample for an array
    hands = datasample(allcards,1);
    
    
    %% Removing the picked card from the deck
    allcards(allcards==hands) = []; 
    allcards = [allcards,hands,hands,hands]; 
    
    
    %% Picking up one more card from the deck
    hands = [hands,datasample(allcards,1)];
    
    
    %% Checking whether hands is a pair
    if hands(1) == hands(2)
        p(k) = 1;
    end
    pp(i) = (sum(p)/i) * 100;
end
plot([1:k],pp,'r.')
percentage = pp(end)
end

Answer 1

[Guillaume]

Monte carlo works fine. Clearly you have a bug.

Why are you shuffling (with randperm) each suit before concatenating them? Would the order of the card in each suit matter when you draw one card at random?

Your allcards array is basically the numbers 1:13 repeated 4 times. You then draw one card (randsample is a bit overkill for this), so a number 1 to 13, and remove all instances of that number hence all 4 cards with the same value from the set, before putting back 3 of them at the end. Finally you draw one card again. Why such complicated logic? Why not draw two cards at once without replacement and be done with it?

A much simpler way to implement that simulation:

p = zeros(1, k);
allcards = repelem(1:13, 4);  %order of cards is irrelevant
for draw = 1:k
    hand = allcards(randperm(52, 2));  %draw two distinct cards at random
    p = diff(hand) == 0;
end
pp = cumsum(p) ./ (1:k);
plot(p*100);

which indeed converges to 100 x 3/51

edit: And I've finally spotted the bug in your code because when I looked at the p that your code generates it's all 0s except for the last element

if hands(1) == hands(2)
    p(k) = 1;

with k constant and equal to the last element of p regardless of the loop index.

[Day Hong Kim]

Thank you so much Guillaume! Due to my lack of knowledge, I only know some built-in functions. I will study more. Have a nice Thanksgiving!