盘比赛赔率仿真（MATLAB）

Question

我最近发现了巨大的卡附带 - SET。简单地说，有81卡的四个特点：符号（椭圆形，波浪线或菱形），颜色（红色，紫色或绿色），号（一个，两个或三个）和阴影（固体，条纹或开放）。其任务是从选定的12张卡片中找出一组3张卡片，其中每张卡片上的四个特征中的每一个都是相同的，或者每张卡片上的所有不同（不是2 + 1组合）。

我编写它在MATLAB找到一个解决方案，并估计其在随机抽取卡片一套的几率。

这里是我的代码来估计赔率：

%% initialization 
K = 12; % cards to draw 
NF = 4; % number of features (usually 3 or 4) 
setallcards = unique(nchoosek(repmat(1:3,1,NF),NF),'rows'); % all cards: rows - cards, columns - features 
setallcomb = nchoosek(1:K,3); % index of all combinations of K cards by 3 

%% test 
tic 
NIter=1e2; % number of test iterations 
setexists = 0; % test results holder 
% C = progress('init'); % if you have progress function from FileExchange 
for d = 1:NIter 
% C = progress(C,d/NIter);  

% cards for current test 
setdrawncardidx = randi(size(setallcards,1),K,1); 
setdrawncards = setallcards(setdrawncardidx,:); 

% find all sets in current test iteration 
for setcombidx = 1:size(setallcomb,1) 
    setcomb = setdrawncards(setallcomb(setcombidx,:),:); 
    if all(arrayfun(@(x) numel(unique(setcomb(:,x))), 1:NF)~=2) % test one combination 
     setexists = setexists + 1; 
     break % to find only the first set 
    end 
end 
end 
fprintf('Set:NoSet = %g:%g = %g:1\n', setexists, NIter-setexists, setexists/(NIter-setexists)) 
toc 

100-1000迭代的速度快，但要小心了。在我的家用电脑上进行一百万次迭代大约需要15个小时。无论如何，有12张牌和4个特征，我有13：1的套牌。这实际上是一个问题。说明书说这个数字应该是33：1。最近由Peter Norvig确认。他提供了Python代码，但我还没有测试它。

所以你能找到一个错误？欢迎任何关于性能改进的意见。

Answer 1

我在查看代码之前解决了编写自己的实现的问题。我第一次尝试是非常相似，你已经有了:)

%# some parameters 
NUM_ITER = 100000; %# number of simulations to run 
DRAW_SZ = 12;  %# number of cards we are dealing 
SET_SZ = 3;   %# number of cards in a set 
FEAT_NUM = 4;  %# number of features (symbol,color,number,shading) 
FEAT_SZ = 3;  %# number of values per feature (eg: red/purple/green, ...) 

%# cards features 
features = { 
    'oval' 'squiggle' 'diamond' ; %# symbol 
    'red' 'purple' 'green' ;   %# color 
    'one' 'two' 'three' ;   %# number 
    'solid' 'striped' 'open'   %# shading 
}; 
fIdx = arrayfun(@(k) grp2idx(features(k,:)), 1:FEAT_NUM, 'UniformOutput',0); 

%# list of all cards. Each card: [symbol,color,number,shading] 
[W X Y Z] = ndgrid(fIdx{:}); 
cards = [W(:) X(:) Y(:) Z(:)]; 

%# all possible sets: choose 3 from 12 
setsInd = nchoosek(1:DRAW_SZ,SET_SZ); 

%# count number of valid sets in random draws of 12 cards 
counterValidSet = 0; 
for i=1:NUM_ITER 
    %# pick 12 cards 
    ord = randperm(size(cards,1)); 
    cardsDrawn = cards(ord(1:DRAW_SZ),:); 

    %# check for valid sets: features are all the same or all different 
    for s=1:size(setsInd,1) 
     %# set of 3 cards 
     set = cardsDrawn(setsInd(s,:),:); 

     %# check if set is valid 
     count = arrayfun(@(k) numel(unique(set(:,k))), 1:FEAT_NUM); 
     isValid = (count==1|count==3); 

     %# increment counter 
     if isValid 
      counterValidSet = counterValidSet + 1; 
      break   %# break early if found valid set among candidates 
     end 
    end 
end 

%# ratio of found-to-notfound 
fprintf('Size=%d, Set=%d, NoSet=%d, Set:NoSet=%g\n', ... 
    DRAW_SZ, counterValidSet, (NUM_ITER-counterValidSet), ... 
    counterValidSet/(NUM_ITER-counterValidSet)) 

使用Profiler来发现热点后，一些改进可以主要由由早期break'ing出循环的可能时。主要的瓶颈是对UNIQUE函数的调用。我们检查有效集合的上述两行可以改写为：

%# check if set is valid 
isValid = true; 
for k=1:FEAT_NUM 
    count = numel(unique(set(:,k))); 
    if count~=1 && count~=3 
     isValid = false; 
     break %# break early if one of the features doesnt meet conditions 
    end 
end 

不幸的是，对于较大的仿真，仿真仍然很慢。因此，我的下一个解决方案就是矢量化版本，每次迭代时，我们都会从12张牌的牌手中构建出所有可能的3张牌的单个矩阵。对于所有这些候选集合，我们使用逻辑向量来指示出现哪些特征，从而避免调用UNIQUE/NUMEL（我们希望每个卡片上的特征都相同或不同）。我承认代码现在不易读，也很难遵循（因此我发布了两个版本进行比较）。原因是我试图尽可能优化代码，因此每个迭代循环都是完全向量化的。下面是最终代码：

%# some parameters 
NUM_ITER = 100000; %# number of simulations to run 
DRAW_SZ = 12;  %# number of cards we are dealing 
SET_SZ = 3;   %# number of cards in a set 
FEAT_NUM = 4;  %# number of features (symbol,color,number,shading) 
FEAT_SZ = 3;  %# number of values per feature (eg: red/purple/green, ...) 

%# cards features 
features = { 
    'oval' 'squiggle' 'diamond' ; %# symbol 
    'red' 'purple' 'green' ;   %# color 
    'one' 'two' 'three' ;   %# number 
    'solid' 'striped' 'open'   %# shading 
}; 
fIdx = arrayfun(@(k) grp2idx(features(k,:)), 1:FEAT_NUM, 'UniformOutput',0); 

%# list of all cards. Each card: [symbol,color,number,shading] 
[W X Y Z] = ndgrid(fIdx{:}); 
cards = [W(:) X(:) Y(:) Z(:)]; 

%# all possible sets: choose 3 from 12 
setsInd = nchoosek(1:DRAW_SZ,SET_SZ); 

%# optimizations: some calculations taken out of the loop 
ss = setsInd(:); 
set_sz2 = numel(ss)*FEAT_NUM/SET_SZ; 
col = repmat(1:set_sz2,SET_SZ,1); 
col = FEAT_SZ.*(col(:)-1); 
M = false(FEAT_SZ,set_sz2); 

%# progress indication 
%#hWait = waitbar(0./NUM_ITER, 'Simulation...'); 

%# count number of valid sets in random draws of 12 cards 
counterValidSet = 0; 
for i=1:NUM_ITER 
    %# update progress 
    %#waitbar(i./NUM_ITER, hWait); 

    %# pick 12 cards 
    ord = randperm(size(cards,1)); 
    cardsDrawn = cards(ord(1:DRAW_SZ),:); 

    %# put all possible sets of 3 cards next to each other 
    set = reshape(cardsDrawn(ss,:)',[],SET_SZ)'; 
    set = set(:); 

    %# check for valid sets: features are all the same or all different 
    M(:) = false;   %# if using PARFOR, it will complain about this 
    M(set+col) = true; 
    isValid = all(reshape(sum(M)~=2,FEAT_NUM,[])); 

    %# increment counter if there is at least one valid set in all candidates 
    if any(isValid) 
     counterValidSet = counterValidSet + 1; 
    end 
end 

%# ratio of found-to-notfound 
fprintf('Size=%d, Set=%d, NoSet=%d, Set:NoSet=%g\n', ... 
    DRAW_SZ, counterValidSet, (NUM_ITER-counterValidSet), ... 
    counterValidSet/(NUM_ITER-counterValidSet)) 

%# close progress bar 
%#close(hWait) 

如果你有并行处理工具箱，你可以很容易地并行PARFOR替代for循环平原（您可能要再次移动矩阵M的初始化循环内：与M = false(FEAT_SZ,set_sz2);替换M(:) = false;）

这里是50000个模拟（PARFOR 2个本地实例的池使用的）的一些示例输出：

» tic, SET_game2, toc 
Size=12, Set=48376, NoSet=1624, Set:NoSet=29.7882 
Elapsed time is 5.653933 seconds. 

» tic, SET_game2, toc 
Size=15, Set=49981, NoSet=19, Set:NoSet=2630.58 
Elapsed time is 9.414917 seconds. 

并与一百万次迭代（PARFOR 12，无PARFOR 15）：

» tic, SET_game2, toc 
Size=12, Set=967516, NoSet=32484, Set:NoSet=29.7844 
Elapsed time is 110.719903 seconds. 

» tic, SET_game2, toc 
Size=15, Set=999630, NoSet=370, Set:NoSet=2701.7 
Elapsed time is 372.110412 seconds. 

的比值比同意由Peter Norvig结果报告如下。

Answer 2

这里有一个量化版本，其中100万手可以在一分钟左右来计算。我得到了大约28：1的比例，所以在找到“所有不同的”组合后可能还会有一些小小的变化。我的猜测是，这也是你的解决方案遇到的问题。

%# initialization 
K = 12; %# cards to draw 
NF = 4; %# number of features (this is hard-coded to 4) 
nIter = 100000; %# number of iterations 

%# each card has four features. This means that a card can be represented 
%# by a coordinate in 4D space. A set is a full row, column, etc in 4D 
%# space. We can even parallelize the iterations, at least as long as we 
%# have RAM (each hand costs 81 bytes) 
%# make card space - one dimension per feature, plus one for the iterations 
cardSpace = false(3,3,3,3,nIter); 

%# To draw cards, we put K trues into each cardSpace. I can't think of a 
%# good, fast way to draw exactly K cards that doesn't involve calling 
%# unique 
for i=1:nIter 
    shuffle = randperm(81) + (i-1) * 81; 
    cardSpace(shuffle(1:K)) = true; 
end 

%# to test, all we have to do is check whether there is any row, column, 
%# with all 1's 
isEqual = squeeze(any(any(any(all(cardSpace,1),2),3),4) | ... 
    any(any(any(all(cardSpace,2),1),3),4) | ... 
    any(any(any(all(cardSpace,3),2),1),4) | ... 
    any(any(any(all(cardSpace,4),2),3),1)); 
%# to get a set of 3 cards where all symbols are different, we require that 
%# no 'sub-volume' is completely empty - there may be something wrong with this 
%# but since my test looked ok, I'm not going to investigate on Friday night 
isDifferent = squeeze(~any(all(all(all(~cardSpace,1),2),3),4) & ... 
    ~any(all(all(all(~cardSpace,1),2),4),3) & ... 
    ~any(all(all(all(~cardSpace,1),3),4),2) & ... 
    ~any(all(all(all(~cardSpace,4),2),3),1)); 

isSet = isEqual | isDifferent; 

%# find the odds 
fprintf('odds are %5.2f:1\n',sum(isSet)/(nIter-sum(isSet))) 

Answer 3

我发现我的错误。感谢Jonas与RANDPERM的提示。

我用RANDI以随机抽取的ķ卡，但有大约50％的机会得到重复，即使在12张牌。当我用randperm替换这一行时，我得到了33.8：1的10000次迭代，非常接近指令手册中的数字。

setdrawncardidx = randperm(81); 
setdrawncardidx = setdrawncardidx(1:K); 

无论如何，看到其他解决问题的方法会很有趣。

Answer 4

我确定我的这些几率的计算有些问题，因为其他几个人已经通过模拟证实了它接近于33：1，正如在说明中那样，但以下逻辑有什么问题？

对于12张随机卡，有三张牌220个可能的组合（12！/（9！3！）= 220）。三张卡片的每个组合都有1/79的机会，所以三张任意卡片不会成为一组的概率为78/79。所以如果你检查了所有的220种组合，并且每个组合都有78/79的机会，那么你没有找到一套检查所有可能组合的机会将被78/79提升到220次方，或0.0606，这大约是。 17：1的赔率。

我必须错过什么......？

Christopher