MATLAB Answers

Monte Carlo Method Fails

4 views (last 30 days)
Day Hong Kim
Day Hong Kim on 26 Nov 2019
Commented: Day Hong Kim on 26 Nov 2019
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;
pp(i) = (sum(p)/i) * 100;
percentage = pp(end)


Sign in to comment.

Accepted Answer

Guillaume on 26 Nov 2019
Edited: Guillaume on 26 Nov 2019
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;
pp = cumsum(p) ./ (1:k);
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.

  1 Comment

Day Hong Kim
Day Hong Kim on 26 Nov 2019
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!

Sign in to comment.

More Answers (0)

Sign in to answer this question.