After 2*(n-1) draws, there will be either two black beans or one white bean remaining. All cases of just one white bean and nothing else remaining occur at this stage. Thus we keep track of the chain of probabilities of each possible state of number of black and white beans, for 2*(n-1) draws and find the probability of state (0,1) at that point.
DefDbl A-Z
Dim crlf$
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()
Form1.Visible = True
Text1.Text = ""
crlf = Chr$(13) + Chr$(10)
For n = 2 To 10
ReDim oldstate(n, n, 1) ' black,white
oldstate(0, n, 0) = 1 ' numerator
oldstate(0, n, 1) = 1 ' denominator
For draw = 1 To 2 * (n - 1)
ReDim newstate(n, n, 1)
For black = 0 To n
For white = 0 To n
If oldstate(black, white, 0) > 0 Then
If white > 0 Then If newstate(black + 1, white - 1, 1) = 0 Then newstate(black + 1, white - 1, 1) = 1
If black > 0 Then If newstate(black - 1, white, 1) = 0 Then newstate(black - 1, white, 1) = 1
If white > 0 Then ' choose white
lcd = lcm((black + white) * oldstate(black, white, 1), newstate(black + 1, white - 1, 1))
num = newstate(black + 1, white - 1, 0) * lcd / newstate(black + 1, white - 1, 1) + white * oldstate(black, white, 0) * lcd / oldstate(black, white, 1) / (black + white)
den = lcd
g = gcd(num, den)
newstate(black + 1, white - 1, 0) = num / g
newstate(black + 1, white - 1, 1) = den / g
End If
If black > 0 Then ' choose black
lcd = lcm((black + white) * oldstate(black, white, 1), newstate(black - 1, white, 1))
num = newstate(black - 1, white, 0) * lcd / newstate(black - 1, white, 1) + black * oldstate(black, white, 0) * lcd / oldstate(black, white, 1) / (black + white)
den = lcd
g = gcd(num, den)
newstate(black - 1, white, 0) = num / g
newstate(black - 1, white, 1) = den / g
End If
End If
Next white
Next black
For black = 0 To n
For white = 0 To n
oldstate(black, white, 0) = newstate(black, white, 0)
oldstate(black, white, 1) = newstate(black, white, 1)
Next
Next
' If n = 5 Then
' For black = 0 To n
' For white = 0 To n
' If newstate(black, white, 1) = 0 Then newstate(black, white, 1) = 1
' Text1.Text = Text1.Text & mform(newstate(black, white, 0) / newstate(black, white, 1), "##0.000000")
' Next
' Text1.Text = Text1.Text & crlf
' Next
' Text1.Text = Text1.Text & crlf
' End If
Next draw
Text1.Text = Text1.Text & n & Str(newstate(0, 1, 0)) & "/" & (newstate(0, 1, 1)) & " ~= " & newstate(0, 1, 0) / newstate(0, 1, 1) & crlf
Next n
Text1.Text = Text1.Text & crlf & " done"
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
Function lcm(a, b)
lcm = a * b / gcd(a, b)
End Function
finding
n probability
2 1/2 ~= 0.5
3 7/18 ~= 0.388888888888889
4 97/288 ~= 0.336805555555556
5 54997/180000 ~= 0.305538888888889
6 460463/1620000 ~= 0.284236419753086
7 51185279267/190591380000 ~= 0.26856030564971
8 200170674968477/780662292480000 ~= 0.256411353406832
9 2.8197972562852E+15/1.143286897536E+16 ~= 0.246639514750182
10 2.75730948120933E+15/1.15582791307704E+16 ~= 0.238557094011412
So for n=5, the solution is 54997/180000 ~= 0.305538888888889.
Up to n=8, the denominators are A121564 (Denominator of Sum[i=1..n] i!/(i^i)):
n a(n)
1 1
2 2
3 18
4 288
5 180000
6 1620000
7 190591380000
8 780662292480000
9 46097327708651520000
10 28810829817907200000000
11 747278726094210559615027200000000
12 747278726094210559615027200000000
The mismatch beyond n=8 might be the result of the limited precision available interfering with GCD calculation; note the exponential default display format.
But a search for the numerators is unsuccessful.
Simulation:
n hits trials
2 49877 100000
3 39145 100000
4 33726 100000
5 30691 100000
6 28302 100000
7 27000 100000
8 25567 100000
9 24486 100000
10 23847 100000
11 22950 100000
12 22562 100000
13 22044 100000
14 21667 100000
15 21159 100000
16 20870 100000
17 20597 100000
18 20247 100000
19 19922 100000
20 19651 100000
...
40 17015 100000
50 15927 100000
...
100 13968 100000
...
1000 10584 100000
...
10000 8102 100000
Indeed the probability seems it must converge: that is, it seems to be monotonically decreasing and there's a least value that it can be, zero. If it indeed is monotonically decreasing it must coverge to some finite value, possibly zero.
Simulation from:
DefDbl A-Z
Dim crlf$
Private Sub Form_Load()
Form1.Visible = True
Text1.Text = ""
crlf = Chr$(13) + Chr$(10)
For n = 2 To 20
Randomize Timer
hitct = 0: ct = 0
For tr = 1 To 100000
white = n: black = 0
Do While white + black > 0
r = Int(Rnd(1) * (white + black) + 1)
If r <= white Then
black = black + 1: white = white - 1
Else
black = black - 1
End If
If white = 1 And black = 0 Then hitct = hitct + 1: Exit Do
Loop
ct = ct + 1
Next tr
Text1.Text = Text1.Text & n & Str(hitct) & Str(ct) & crlf
DoEvents
Next n
Text1.Text = Text1.Text & crlf & " done"
End Sub
|
|
Posted by Charlie
on 2015-07-13 18:49:13 |
computer aided solution
