A standard six-sided die is to be rolled repeatedly until a side appears a number of times equal to its number. In other words until the n-th n appears.
Let P(n)=the probability the game terminates with the n-th n.
Find the distribution of n.
Feel free to generalize for m sides.
Warning: I have not managed this past m=4.
n prob(decimal) prob (exact)
1 0.6881516453401 8987738063/13060694016
2 0.2125194475337 925217159/4353564672
3 0.0685578584217 895413211/13060694016
4 0.0218963008044 285980885/13060694016
5 0.0068166367662 9892223/1451188224
6 0.0020581121476 26880373/13060694016
1.0000000010136
144.252
The UBASIC program
10 a= 8987738063//13060694016
20 b= 925217159//4353564672
30 c= 895413211//13060694016
40 d= 285980885//13060694016
50 e= 9892223//1451188224
60 f= 26880373//13060694016
70 print a+b+c+d+e+f
verifies the total of the exact fractions is exactly 1. The decimals, as can be seen from the decimal total, cannot be trusted beyond 8 places, as they came from many small pieces; in fact 572,493,416 additions, as determined by adding a counter to the program, of the probabilities went into the calculation: one for each possible sequence leading to a positive result. The 144.252 is the number of seconds the calculation took (just counting the number of additions added 40 seconds to this total).
The calculation was in the form of a recursive routine, that fortunately doesn't have to go through all 6^16 possibilities (just the above-mentioned over half a billion possibilities) of the first 16 throws (at which point you're assured that some number will hit, as at 15 throws, you could have no 1's, one 2, two 3's, three 4's, four 5's and five 6's). As soon as you account for rolling a 1, you need not go down the path of whatever may follow; similarly for the second 2, etc.:
DefDbl A-Z
Dim crlf$, currprob, currnum, p4no(6), had(6), hist(16), num(6)
Function mform$(x, t$)
a$ = Format$(x, t$)
If Len(a$) < Len(t$) Then a$ = Space$(Len(t$) - Len(a$)) & a$
mform$ = a$
End Function
Private Sub Form_Load()
t = Timer
Text1.Text = ""
crlf$ = Chr(13) + Chr(10)
Form1.Visible = True
currprob = 1
currnum = Int(6 ^ 16 + 0.5) ' first numerator, to make fraction = 1
doNext 1
For i = 1 To 6
Text1.Text = Text1.Text & i & mform(p4no(i), " 0.0000000000000")
g = gcd(num(i), Int(6 ^ 16 + 0.5))
Text1.Text = Text1.Text & " " & num(i) / g & "/" & Int(6 ^ 16 + 0.5) / g & crlf
totl = totl + p4no(i)
Next
Text1.Text = Text1.Text & mform(totl, "0.0000000000000") & crlf
Text1.Text = Text1.Text & Timer - t
End Sub
Sub doNext(wh)
prevprob = currprob
currprob = currprob / 6
currnum = Int(currnum / 6 + 0.5)
For choice = 1 To 6
had(choice) = had(choice) + 1
hist(wh) = choice
If wh = 2 Then Text2.Text = hist(1) & Str(hist(2)): DoEvents
If had(choice) = choice Then
p4no(choice) = p4no(choice) + currprob
num(choice) = num(choice) + currnum
Else
doNext wh + 1
End If
had(choice) = had(choice) - 1
Next
currprob = prevprob
currnum = currnum * 6
End Sub
Function gcd(a, b)
x = a: y = b
Do
q = Int(x / y)
z = x - q * y
x = y: y = z
Loop Until z = 0
gcd = x
End Function
Simulation verifies
1 6881809 .6881809
2 2123866 .2123866
3 686921 .0686921
4 218414 .0218414
5 68281 .0068281
6 20709 .0020709
where the code
Randomize Timer
For trial = 1 To 1000000
ReDim h(6)
Do
r = Int(6 * Rnd(1) + 1)
h(r) = h(r) + 1
If h(r) = r Then
ct(r) = ct(r) + 1
trct = trct + 1
Exit Do
End If
Loop
Next
For r = 1 To 6
Text1.Text = Text1.Text & r & Str(ct(r)) & Str(ct(r) / trct) & crlf
Next
was executed 10 times (in response to manually clicking a button each time) with the counts accumulating to simulate 10 million trials.
|
|
Posted by Charlie
on 2015-07-01 17:00:33 |
computer aided solution
