[VA] SRC #012a - Then and Now: Probability

[Valentin Albillo]

  
Hi, all,

Thank you very much for your interest in my SRC #012a, and in particular thanks to Vincent Weber, Fernando del Rey, J-F Garnier, C.Ret, Werner, Albert Chan, PeterP, ijabbot, ttw, rprosperi, rawi, and pier4r for your solutions and/or comments, much appreciated and very fine efforts throughout.

In my OP I mentioned "takes random steps" with the idea of misleading people into thinking that Monte Carlo simulations would be the way to go, but as some of you discovered upon attempting the feat, even a large numer of simulations would get just one correct digit for the toy case, let alone the (30, 60) case, as demonstrated by J-F Garnier's awesome little piece of code.

Getting the 10-12 digits asked for required a fully deterministic method and there are at least two ways to proceed: (a) for each point, identify its neighbors and take the due part from them, or (b) for each location, identify its neighbors and distribute the point's probability among them. It's a zero-sum game and the total probability is the same at the beginning (P=1 for the top location, 0 for all the rest), after every step and thus at the end.

I tried both approaches (a) and (b) above and both worked alright but approach (b), distributing the probability among the neighbors, proved to be faster so that's what my original solution uses. Though tempting and perfectly adequate for the (30, 60) case as originally stated, I don't use symmetry for the reasons given in the fourth Note below, mainly because it would lose generality, which I wanted to preserve so that I could also solve arbitrary starting positions (non-symmetric, one or more as long as their initial peobabilities added up to exactly 1).

So, this is my original solution, a 13-line, 683-byte program. For convenience, it uses two matrix keywords from the Math ROM, namely MAT..ZER and MAT=, easily replaced by trivial loops if the Math ROM's not available):

PROBLEM1 (683 bytes)

10 DESTROY ALL @ OPTION BASE 1 @ INPUT "Rows,Steps=";M,N @ DIM A(M,M),B(M,M) @ SETTIME 0
15 A(1,1)=1 @ W=M-1 @ FOR I=1 TO N @ MAT B=ZER
20 P=A(1,1)/2 @ IF P THEN B(2,1)=B(2,1)+P @ B(2,2)=B(2,2)+P
25 P=A(M,1)/2 @ IF P THEN B(W,1)=B(W,1)+P @ B(M,2)=B(M,2)+P
30 P=A(M,M)/2 @ IF P THEN B(W,W)=B(W,W)+P @ B(M,W)=B(M,W)+P
35 FOR X=2 TO W @ U=X-1 @ V=X+1 @ P=A(X,1)/4 @ Q=A(X,X)/4 @ R=A(M,X)/4
40 IF P THEN B(U,1)=B(U,1)+P @ B(V,1)=B(V,1)+P @ B(X,2)=B(X,2)+P @ B(V,2)=B(V,2)+P
45 IF Q THEN B(U,U)=B(U,U)+Q @ B(V,V)=B(V,V)+Q @ B(X,U)=B(X,U)+Q @ B(V,X)=B(V,X)+Q
50 IF R THEN B(M,U)=B(M,U)+R @ B(M,V)=B(M,V)+R @ B(W,U)=B(W,U)+R @ B(W,X)=B(W,X)+R
55 NEXT X @ FOR X=3 TO W @ U=X-1 @ V=X+1 @ FOR Y=2 TO U @ R=Y-1 @ S=Y+1 @ P=A(X,Y)/6
60 IF P THEN B(U,R)=B(U,R)+P @ B(U,Y)=B(U,Y)+P @ B(X,R)=B(X,R)+P
65 IF P THEN B(X,S)=B(X,S)+P @ B(V,Y)=B(V,Y)+P @ B(V,S)=B(V,S)+P
70 NEXT Y @ NEXT X @ MAT A=B @ NEXT I @ P=0 @ FOR Y=1 TO M @ P=P+A(M,Y) @ NEXT Y @ DISP P;TIME

>STD @ RUN
Rows,Steps= 5,4 [ENDLINE] -> 7.98611111112E-2 (=23/288) .03"

After running the program, we can display the just-built probability matrix (which contains the probability that the man is in each of the 5*(5+1)/2 = 15 points in the grid after he's walked the 4 steps in this toy case, which is small enough that we can get it in exact rational form (requires the JPC ROM) right from the command line:

>FOR I=1 TO 5 @ FOR J=1 TO I @ FRAC$(A(I,J));" "; @ NEXT J @@ NEXT I


                         11/96
                   49/384     49/384
           113/1152     101/576     113/1152
      35/1152       17/288     17/288       35/1152
1/128       23/1152      7/288      23/1152       1/128

Now for the big (30, 60) case:

>RUN
Rows,Steps= 30,60 [ENDLINE] -> 9.51234350207E-6   27.46"

We can now use the probability matrix right from the command line to answer all kinds of questions, e.g.:

(a) check the matrix correctness by adding up the probabilities for all points:

>S=0 @ FOR I=1 TO M @ FOR J=1 TO I @ S=S+A(I,J) @ NEXT J @ NEXT I @ S

      1.00000000002

which gives 1 as it should, after all the man must be in one of the 465 points.

(b) probability that he ends in the left border:

>L=0 @ FOR I=1 TO M @ L=L+A(I,1) @ NEXT I @ L

      .117372130231

which is 12338.9289091 times more probable than ending in the bottom row.

(c) probability that he ends in any border as compared to ending inside the grid:

>D=0 @ FOR I=1 TO M @ D=D+A(M,I) @ NEXT I @ D

      9.51234350207E-6 (down edge = bottom row)

>R=L @ L+R+D

      .234753772806

so he ends in a border 23.48% of the time and thus 76.52% of the time he ends inside the grid, so he's 3.26 times more likely to end inside the grid than in a border.

(d) probability that he ends in each of the corners or in any of them:

>A(1,1);A(M,1);A(M,M) @ A(1,1)+A(M,1)+A(M,M)

      8.30321536116E-3   2.21000524869E-8   2.21000524869E-8 (symmetric)

      8.30325956126E-3

(e) probability that he ends up in any of the first N rows:

>FOR N=1 TO M @ P=0 @ FOR I=1 TO N @ FOR J=1 TO I @ P=P+A(I,J) @ NEXT J @ NEXT I @ N;P @ NEXT N


 N      P
---------------------------------------------
 1   0.00830      top corner
 2   0.04100      top two rows
       ...
 5   0.25643      top 5 rows
       ...
 7   0.44905      top 7 rows, less than 50%
 8   0.54414      top 8 rows, more than 50%
       ...
10   0.71182      top 10 rows, 71%
15   0.94283      top 15 rows, 94%
20   0.99430      top 20 rows, 99%
25   0.99972      top 25 rows, almost sure
29   0.99999      top 29 rows, 99.999% certain
30   1.00000      all 30 rows, 100% certain.

Notes:

• All runtimes are for a virtual HP-71B (go71b) running on a mid-range Samsung Galaxy Tab A 6 tablet (Android). A physical HP-71B should take ~58 min for the (30, 60) case.
    
• The (30, 60) case needs 14,533 bytes of RAM to run. If the matrices are dimensioned as SHORT then it requires only 8,245 bytes of RAM but it runs slower and gives only 5 correct digits save 2 ulp (P=0.0000095121), instead of 12 correct digits.
    
• The program checks whether a location's probability is currently zero, in which case it skips it when distributing the probability among its neighbors. This saves about 20% running time for the (30,60) case.
    
• Another pretty obvious optimization would be to take advantage of the symmetry, as the man's starting location is at the top corner; this would cut running time and required RAM roughly in half but it loses generality, i.e.: arbitrary non-symmetrically placed starting positions (one or more) would not run at all, and it would complicate the coding somewhat; also, for a one-time program which already runs suitably fast it's really a moot point.
    
• It's a real pity that the HP-71B BASIC dialect doesn't include the += , etc., operators, as it would make my program that much shorter and faster, i.e. line
  
  
  50 IF R THEN B(M,U)=B(M,U)+R @ B(M,V)=B(M,V)+R @ B(W,U)=B(W,U)+R @ B(W,X)=B(W,X)+R
  
  would become
  
  
  50 IF R THEN B(M,U)+=R @ B(M,V)+=R @ B(W,U)+=R @ B(W,X)+=R
  
  and so on. Perhaps J-F Garnier could do something about it.
    
• The exact value is (465 points updated per step * 60 steps = ~ 200,000 arithmetic operations in all):
  
  
     3722200777884626618385530906788866022689096963173522895529  
  391302183008102676318141068027642364938466415279908207464546304
  
  = 9.5123435020743320973531347375383316350767310071737e−6 (50 digits)

Frankly, I thought that some accomplished RPL programmers would use the multiprecision capabilities of their advanced RPL models to compute the exact rational solution above but no such luck ...


Thanks to all of you for your interest, solutions and comments. A number of you produced correct solutions for this Problem 1 but it's the easiest of the lot and I wonder how many of you will achieve a perfect 6 for 6 score by the time all six parts of SRC #012 are over ...

See you in Problem 2 featuring soon.
V.