Prob of not drawing card in x trials...help

[R Gibert]

In my Monte Carlo program P(at least 1 not drawn) is 100% for 90 shift/trials. If I increase the number of shift/trials to 120, I get 99.789%. With 180 shift/trials, I get 81.699%. Here is my Python program:
<font class="small">Code:</font><hr /><pre>
from random import *

n = 100000
t = 90
ALL = 2**52 - 1
c = 0
for j in range(n):
m = 0
for i in range(t):
d = randint(0,51)
m |= (1 &lt;&lt; d)
if m == ALL:
c += 1
print 1 - float(c)/n
</pre><hr />
I am pretty tired, so I must have screwed up somewhere.

However, speaking in general, I trust these simple Monte Carlo programs more than tricky calcs. I realize that for some pathological cases, Monte Carlo programs are not to be trusted, but this case shouldn't be one of them.

One possibility that just occurs to me is that your calc got counterfeited by overflow. I haven't gone over your calc in detail to say for sure that this is what is going on, but you might look into it.

My approach sidesteps this, but in any case this is seldom a problem in Python as it automatically handles huge numbers gracefully e.g 1000! is calculated by:

reduce(lambda x,y:x*y, range(1,1000))
40238726007709377354370243392300398571937486421071 46325
43799910429938512398629020592044208486969404800479 98861
01971960586316668729948085589013238296699445909974 24504
08707375991882362772718873251977950595099527612087 49754
62497043601418278094646496291056393887437886487337 11918
10458257836478499770124766328898359557354325131853 23958
46307555740911426241747434934755342864657661166779 73966
68820291207379143853719588249808126867838374559731 74613
60853795345242215865932019280908782973084313928444 03281
23155861103697680135730421616874760967587134831202 54785
89320767169132448426236131412508780208000261683151 02734
18279777047846358681701643650241536913982812648102 13092
76124489635992870511496497541990934222156683257208 08213
33186116811553615836546984046708975602900950537616 47584
77284218896796462449451607653534081989013854424879 84959
95331910172335555660213945039973628075013783761530 71277
61926849034352625200015888535147331611702103968175 92151
09077880193931781141945452572238655414610628921879 60223
83897147608850627686296714667469756291123408243920 81601
53780889893964518263243671616762179168909779911903 75403
12746222899880051954444142820121873617459926429565 81746
62830295557029902432415318161721046583203678690611 72601
58783520751516284225540265170483304226143974286933 06169
08979684825901254583271682264580665267699586526822 72807
07578139185817888965220816434834482599326604336766 01769
99612831860788386150279465955131156552036093988180 61213
85586003014356945272242063446317974605946825731037 90084
02443243846565724501440282188525247093519062092902 31364
93273497565513958720559654228749774011413346962715 42284
58623773875382304838656889764619273838149001407673 10446
64025989949022222176590433990188601856652648506179 97023
56193897017860040811889729918311021171229845901641 92106
88843871218556461249607987229085192968193723886426 14839
65738229112312502418664935314397013742853192664987 53372
18940694281434118520158014123344828015051399694290 15348
30776445690990731524332782882698646027898643211390 83506
21709500259738986355427719674282224875758676575234 42202
07573630569498825087968928162753848863396909959826 28095
61214509948717012445164612603790293091208890869420 28510
64018215439945715680594187274899809425474217358240 10636
77404595741785160829230135358081840096996372524230 56085
59037006242712434169090041536901059339838357779394 10970
02775347200000000000000000000000000000000000000000 00000
00000000000000000000000000000000000000000000000000 00000
00000000000000000000000000000000000000000000000000 00000
00000000000000000000000000000000000000000000000000 00000
00000000000000000000000000000000000L

BTW, a common source of overflow is if you calculate C(n, r) from it's most common definition.

[R Gibert]

Okay, I did goof. I misinterpreted the problem.

The good news is that your calc appears correct. Here is my amended program:
<font class="small">Code:</font><hr /><pre>
from random import *

n = 100000
t = 90
k = 52
ALL = 2**k - 1
c = 0
deck = range(52)
for j in range(n):
m = 0
for i in range(t):
d = sample(deck,2)
m |= (1 &lt;&lt; d[0])
m |= (1 &lt;&lt; d[1])
if m == ALL:
c += 1
print 1 - float(c)/n

</pre><hr />
Running this for 90 shift/trials yields: 80.316%, which is in rough agreement.

[R Gibert]

46.311% confirms your other calc with the following more general program:
<font class="small">Code:</font><hr /><pre>
from random import *

n = 100000
t = 90
k = 52
c = 0
deck = range(52)
for j in range(n):
m = [False]*52
for i in range(t):
d = sample(deck,2)
m[d[0]] = True
m[d[1]] = True
if len(filter(lambda x: x == False, m)) &lt; 2:
c += 1
print 1 - float(c)/n
</pre><hr />
I traded efficiency for something more readable. Retaining the old approach was workable, but the bit twiddling would have produced an ugly result.