The following solution was supplied by Marijke van Gans
Second attempt, after Brian Medley spotted an egregious error in my first one.
Let P(m) be the probability that after m moves we're back at the origin (as
a function of that m, but also depending on our choice of epsilon).
It is immediately obvious that the chance of returning after an odd number of
moves m is zero. The flea is on the wrong color square of the chessboard. Let
us for even m set n = m/2.
The moves can be divided in "horizontal" moves (east and west) and "vertical"
ones (north and south), both flavors have probability 1/2 at each move. The
probability that the flea has executed h horizontal moves and m-h vertical
ones after m moves is therefore the binomial
(m C h) 0.5^m
where the notation (m C h), "m choose h" stands for m! / (h! (m-h)!) as usual.
Again, for the flea to return to coordinates (0,0), both coordinates must
be zero simultaneously, so the horizontal moves cancel each other, and the
vertical moves separately also cancel. So we find that not only m must be
even, but also h and m-h separately must be even. Let us set i = h/2.
As far as the vertical moves are concerned, if there were m-h of them, the
probability they all canceled out is
- for odd m-h, 0
- for even m-h (which is 2n-2i), (2n-2i C n-i) 0.5^(2n-2i)
because we must have half of them cancelling the other half.
As far as the horizontal moves are concerned, if there were h of them, the
probability they all canceled out is
- for odd h, 0
- for even h (h = 2i), (2i C i) (0.5+epsilon)^i (0.5-epsilon)^i
P(2n) is the sum over all cases (h all even numbers from 0 to 2n, in other
words i all integers from 0 to n) where each case is the probability that
both the horz. moves cancel and the vert. moves cancel, and each i case
weighted by the prob. that the case occurs. That is, P(2n) =
Sum[i=0...n] { "weight" Prob(2i of the 2n moves were horizontal) *
Prob(2n-2i vert moves cancel) *
Prob(2i horz moves cancel) } =
Sum[i=0...n] { (2n C 2i) 0.5^(2n) *
(2n-2i C n-i) 0.5^(2n-2i) *
(2i C i) (0.5+epsilon)^i (0.5-epsilon)^i }
Note the "weights" don't add to 1 as they exclude all cases where an odd
number out of 2n moves are horizontal, we could include them in our sum
but these terms would then have to have a factor zero anyway.
Gathering factors, P(2n) =
Sum[i=0...n] (2n C 2i) (2n-2i C n-i) (2i C i) 0.5^(4n)
* (1+2epsilon)^i (1-2epsilon)^i
where the 1 +/- 2epsilon arise from dividing 0.5 +/- epsilon by 0.5 thereby
bumping up the count of factors 0.5 to a fixed 4n. In fact, we can take those
factors (1+2epsilon) (1-2epsilon) together as (1-4epsilon^2) which i will call
eta. So we have P(2n) as
Sum[i=0...n] (2n C 2i) (2n-2i C n-i) (2i C i) 0.5^(4n) * eta^i
where eta = (1+2epsilon) (1-2epsilon) = 1 - 4 epsilon^2.
Let's see what we can do with the binomials, or rather, quadrinomials.
(2n)! (2n-2i)! (2i)!
-------------- * ------------- * ----- * 0.5^(4n) * eta^i
(2i)! (2n-2i)! (n-i)! (n-i)! i! i!
simplifies to
(2n)!
------------------- * 0.5^(4n) * eta^i
(n-i)! (n-i)! i! i!
and to recap, P(2n) is the Sum[0 <= i <= n] of that.
Marijke, > (2r)! > ------------------- * 0.5^(4r) * eta^p > (r-p)! (r-p)! p! p! Nice derivation. >Now E is of the nature of a simple product (of the probabilities that the >flea doesn't return at the first step AND doesn't return at the second >step AND not at the third step AND...) But are these probabilities independent? Intuitively, the flea has a better chance of being at the start in, say, two time-steps from now _if_ we know it's there now than if we don't know that. I used your program to check the probabilities of return after some small numbers of steps and they were consistent with outputs from my program, so I'm sure your formula is correct. Assuming I'm right in thinking the probs aren't independent, I'm sure your formula still leads pretty directly to a solution. Introducing new notation, let p(r) = prob of return after exactly 2r steps = your formula q(r) = prob of return _for_the_first_time_ after 2r steps Then q(1)=p(1) and p(r)=sum(from i=1 to r)[q(i)p(r-i)] which will allow q(r) to be found as long as the previous q(i)'s have been stored. Then add up the q's.
Brian,
> But are these probabilities independent? Intuitively, the flea has a
better chance of being at the start in, say, two time-steps from now _if_
we know it's there now than if we don't know that. <
Aaaaargh! <rips out plucks of hair.... then decides to read on first>
> which will allow q(r) to be found as long as the previous q(i)'s have been
stored.
Then add up the q's. <
Ah! Will work on it some more tonight (Tue).
Duly reassured, i took it up again from there:
So far so good. At this stage, i cannot treat returning to the origin after
2m moves as an independent event from returning after 2m' moves, especially
if m is near m'. Rather, following Brian's suggestion, introduce Q
Q(2n) := Prob(flea returns to origin after 2n moves AND this
is the first time it does so since the time n=0)
It is clear that
Q(0) is not well defined
Q(2) == P(2)
Q(2n) <= P(2n)
For 2n > 0 we have the following relation, again with thanks to Brian: Let 2n
be again a number of moves, and if the flea returns after 2n moves let 2k be
the most recent earlier move when the flea returned. This is always meaningful
for 2n > 0 because if there were no "proper" previous times, we can still use
2k = 0.
Summing over cases (2k in the range 0 <= 2k < 2n) we see that in each case
the term is the product of the probability that we actually got to the origin
at move 2k -- which is P(2k) -- times the probability that we both got to the
origin again 2n-2k moves later AND that there were no intervening times when
we did so -- which is Q(2n-2k), because random walks are without memory, once
we're back at the origin it's just like as if m=0 again.
P(2n > 0) = Sum[0 <= k < n] P(2k) * Q(2n-2k)
The P and Q depend recursively on eachother. We have
P(0) = 1
and by the cases formula
P(2) = P(0) Q(2)
P(4) = P(0) Q(4) + P(2) Q(2)
P(6) = P(0) Q(6) + P(2) Q(4) + P(4) Q(2)
P(8) = P(0) Q(8) + P(2) Q(6) + P(4) Q(4) + P(6) Q(2)
...
from which we can distill the recursion
Q(2) = P(2)
Q(4) = P(4) - P(2) Q(2)
Q(6) = P(6) - P(4) Q(2) - P(2) Q(4)
Q(8) = P(8) - P(6) Q(2) - P(4) Q(4) - P(2) Q(6)
... = ... - ... ... - ... ... - ... ... - ...
The nice thing about the Q's (seen from n=0) is that they describe mutually
excluding cases. If the flea ever returns then the first time it does so is
either at move 2 or 4 or 6 or 8 or... so it's really the sum of the Q's we
are after. We want to rig epsilon (or eta) such that this sum is 1/2.
Calling the sum of all Q's Q, and the sum of all P's except P(0) (if the sum
exists) P, the recursion above gives us (summing all columns)
Q = P - P * Q(2) - P * Q(4) - P * Q(6) - P * Q(8) ...
which is (summing the row)
Q = P - P * Q
Q + PQ = P
Q = P / (1+P)
or, equivalently,
P - PQ = Q
P = Q / (1-Q)
In particular, if Q = 0.5, we must have P = 1.
This observation saves us from actually having to do the recursion each time
to find the Q's from the P's in order to sum the Q's. We can just sum the P's
and take that epsilon for which the sum P approaches 1.
The following UBasic proglet has a loop to interactively home in on eta.
Epsilon, Eta, N, I as in text above.
QNI: quadrinomial coefficient (2n)! / ( [n-i]! [n-i]! i! i! )
calculated for each increasing i by twice /= i, twice *=[n - old i]
QN0: QNI for i = 0 i.e. (2n)! / (n! n!)
calculated from that of previous n by *= (2n-1), twice *= 2n, /= n
EI: eta^i calculated for next i by *= eta
TN: 2^4n calculated for next n by *= 16
PN: Sum[i=0..n] QNI * EI / TN
(calculated by dividing huge sum by huge TN, this way works best for
UBasic with fixed-point accuracy)
P: Sum[n=1...] PN
Cut: lowest value PN we bother with
10 Cut=10^(-15)
100 input "Eta";EtaIn
110 if EtaIn=0 goto 900 ' endgame
120 Eta=EtaIn
200 N=0 ' init n
210 QN0=1 ' init quadrinomial coeff (2*0)!/(0!0!)
220 TN=1 ' init 2^4N
230 P=0
300 inc N
310 QN0*=N+N-1:QN0*=N+N:QN0/=N:QN0/=N
320 TN*=16
330 PN=0
340 QNI=QN0:EI=1:PN=QN0 ' = QNI * EI
350 I=0:J=N
400 inc I
410 QNI*=J
420 QNI/=I
430 QNI*=J
440 QNI/=I
450 dec J
460 EI*=Eta ' eta^I
470 PN+=QNI*EI ' (quadrinomial) * eta^i
480 if J goto 400
500 PN/=TN
510 P+=PN
530 if PN>Cut goto 300
600 print N,P
610 goto 100
900 print "Epsilon =";sqrt(1-Eta)/4
Some preliminary runs show that a cutoff of 10^15 (first example) is as good
as 10^-18 (second example) to within 12 decimal places, and that bumping down
"point" from 8 to 6 or even the default 4 (i.e. 37 to 28 to 19 decimal places
to the right of the decimal point, internally) leaves us enough precision.
run
Words for fractionals 8 (Decimals for display 38)
Eta? .96
1292 1.13380032755507731432316310285080814819
Eta? 0
Epsilon = 0.05
OK
run
Words for fractionals 8 (Decimals for display 38)
Eta? .96
1621 1.13380032755512396179185046002165219978
Eta? 0
Epsilon = 0.05
OK
run
Words for fractionals 6 (Decimals for display 28)
Eta? .96
1621 1.1338003275551239617918504496
Eta? 0
Epsilon = 0.05
OK
run
Words for fractionals 4 (Decimals for display 19)
Eta? .96
1620 1.1338003275551239159
Eta? 0
Epsilon = 0.05
OK
Here's a transcript of the human-computer interaction at point 4, cut 10^-15
run
Eta? .96
1292 1.1338003275550772781
Eta? .95
1038 1.063869184276825339
Eta? .94
868 1.0068636660576407334
Eta? .939
854 1.0017019150950364481
Eta? .9388
851 1.0006798555526157647
Eta? .9387
850 1.0001700938332077058
Eta? .93868
849 1.0000682425793312957
Eta? .93867
849 1.0000173295707334315
Eta? .938667
849 1.0000020573079481553
Eta? .9386667
849 1.0000005301232891355
Eta? .9386666
849 1.0000000210634176628
Eta? .938666596
849 1.0000000007010402924
Eta? .9386665959
849 1.0000000001919808753
Eta? .93866659587
849 1.0000000000392630507
Eta? .93866659586
849 0.9999999999883571081
Eta? .938666595862
849 0.999999999998538296
Eta? .9386665958623
849 1.0000000000000654751
Eta? .93866659586229
849 1.0000000000000145689
Eta? 0
Epsilon = 0.0619139544739865166
OK
Let's throw a few more digits precision at it to check we have the number of
digits we think we have, and just for the heck of it:
point 6
Words for fractionals 6 (Decimals for display 28)
OK
list 10
10 Cut=10^(-25)
run
Eta? .93866659586229
1563 1.0000000000000446415382673716
Eta? .93866659586228
1563 0.9999999999999937355966369573
Eta? .938666595862282
1563 1.0000000000000039167849630396
Eta? .938666595862281
1563 0.9999999999999988261907999984
Eta? .9386665958622812
1563 0.9999999999999998443096326064
Eta? .93866659586228122
1563 0.9999999999999999461215158673
Eta? .93866659586228123
1563 0.9999999999999999970274574978
Eta? .938666595862281231
1563 1.0000000000000000021180516608
Eta? .9386665958622812306
1563 1.0000000000000000000818139957
Eta? 0
Epsilon = 0.0619139544739909428400630502
OK
Checking one last time with yet a bit more wellie, and much higher internal
precision:
point 10
Words for fractionals 10 (Decimals for display 48)
OK
list 10
10 Cut=10^(-30)
OK
run
Eta? .9386665958622812306
1923 1.000000000000000000081817053620054845406966873213
Eta? 0
Epsilon = 0.061913954473990942840063050290634662106957964569
OK
Yes, that looks fine. Finally, check how sensitive epsilon is on eta:
list 900
900 print "Epsilon =";sqrt(1-Eta)/4
OK
Eta = .9386665958622812306000001
OK
goto 900
Epsilon = 0.061913954473990942840062999817360477807160520543
OK
Epsilon = 0.06191395447399094 looks safe as houses to me.
--marijke