 All about flooble | fun stuff | Get a free chatterbox | Free JavaScript | Avatars  perplexus dot info  Bean coloring (Posted on 2015-07-13) Start with a bag containing 5 white beans. Randomly draw one at a time employing the following rule:

If the bean is white, color it black and put it back in the bag;
If the bean is black, keep it out.

What is the probability that at some point there will be a single white bean in the bag?

Does the probability converge, and if so, to what value?

 No Solution Yet Submitted by Jer Rating: 4.0000 (1 votes) Comments: ( Back to comment list | You must be logged in to post comments.) computer aided solution | Comment 1 of 11
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

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   probability2 1/2 ~= 0.53 7/18 ~= 0.3888888888888894 97/288 ~= 0.3368055555555565 54997/180000 ~= 0.3055388888888896 460463/1620000 ~= 0.2842364197530867 51185279267/190591380000 ~= 0.268560305649718 200170674968477/780662292480000 ~= 0.2564113534068329 2.8197972562852E+15/1.143286897536E+16 ~= 0.24663951475018210 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\$

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 Please log in:

 Search: Search body:
Forums (4)