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

[Valentin Albillo]

  
Hi, all,

After 4 weeks to the day, it seems this SRC #012a has run its course, so time for a few additional comments and a few new results.

The aditional comments

Gjermund Skailand Wrote:This has been a very interesting thread. This is a sys-rpl version for HP50g of CReth SRC12a. On an actual HP50g the calculation time for the 30 60 problem is 5 min 21sec. p= 9.51234350207E-6

Thank you for your appreciation. Your SysRPL version looks amazing, kinda assembly language, producing the correct 12-digit result at least 10x faster than a physical HP-71B, which is truly awesome.

Can someone please confirm that the listing is correct and will produce the stated result in the stated time ?

Albert Chan Wrote:For S = R-1, we can treat triangle as without bottom edge (and the 2 corners) [...] It simplified to:

P(R, S=R-1) = 3^(3-R) - 2^(5-2*R)

Very nice exact symbolic result for that particular case, Albert Chan, congratulations !

Normally this could be construed as going against my stated rules but as you previously posted tons of actual HP-71B code, I'm not complaining.

Werner Wrote:Well, to make Valentin's wishes come true, here's an entry that will solve the 30/60 problem on a real 42S, the only vintage RPN calculator able to do it. [...] Estimate of real 42S running time: 3h05m.

Thank you very much, Werner, for taking my wishes into consideration and producing such a fine HP-42S solution. Your running time estimation on a physical HP-42S seems to be 3x slower than my solution running on a physical HP-71B but, as you say, optimization would probably reduce the timing considerably and some compromises had to be made to fit it into the available RAM.

Nevertheless, running your program in Free42 on my Samsung tablet takes just 4.5 seconds, while still within the rules.

The new results

As I've stated oftentimes, one of my main goals is to get people who like vintage HP calculators to not consider them as obsolete gadgets only fit for collecting or nostalgia, with no real place in the real world, but as still useful devices which can indeed be used to solve modern problems and best of all, to improve one's sleuthing and programming abilities while attempting the solution, in the way of "Experimental Mathematics" (EM for short); quoting from Wikipedia:

"Experimental mathematics is an approach to mathematics in which computation is used to investigate mathematical objects and identify properties and patterns."

In what follows, I'll describe my own EM approach to this Problem 1, always using my program (as listed in Post #50) to do the sleuthing. First of all we get this assorted data, in SCI 6 for easier typing:

Rows Steps  Probability     Rows Steps  Probability
---------------------------------------------------
5    4      7.986111e-2     20   20     5.006369e-8
     40     2.666558e-1          200    5.204890e-2
     50     2.666660e-1          1000   6.666566e-2
     60     2.666666e-1          2000   6.666667e-2
     70     2.666667e-1          

10   10     1.317953e-3     30   30     1.289121e-12 
     100    1.317596e-1          60     9.512344e-6
     1000   1.333333e-1          120    1.694782e-3
                                 240    1.531651e-2
                                 300    2.225664e-2
                                 480    3.539453e-2
                                 960    4.368177e-2

and it clearly seems that there's a limit for the value of the probability P as the number of steps increases, which is P_lim(5) = 2.666667e-1 = 4/15,   P_lim(10) = 2/15,   P_lim(20) = 1/15 and though not so fast  P_lim(30) seems to be converging to 4.444444e-2 = 2/45 , so recognizing the obvious pattern we might then conjecture that in the limit we have, for N rows:

                P_lim(N) = 4/(3*N)

Now we can check additional cases (say N = 7, 13, 22, ... rows) to see if the conjectured formula holds, and if it does we can attempt to find a symbolical proof for it, A. Chan-style !

So far this applies to the probability of being in the bottom row at the end of the walk, but what about the probability in the limit of being in a particular, single location as the number of steps grows indefinitely ? Adding this line to my program will display the resulting probability matrix which gives the probability for each and every grid point, using the FRAC$ keyword from the JPC ROM to output exact rational results:

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

Running the program for 5 rows and a sufficiently large number of steps (S=100), we get in SCI 6:

>RUN
       Rows,Steps=5,100 -> 2.666667e-1


            1/30
         1/15 1/15
      1/15  1/10  1/15
   1/15  1/10  1/10  1/15
1/30  1/15  1/15  1/15  1/30

and we see that P(corners) = 1/30, P(edges) = 1/15 = 2/30 and P(inner) = 1/10 = 3/30, so we conjecture that the ratios are

P(corners) : P(edges) : P(inner) = 1 : 2 : 3

which A. Chan also discovered and posted here. As a check, if we run my program for 10 rows, we'll get P(corners), P(edges), P(inner) = 1/135, 2/135, 1/45 = 3/135, further confirming the 1 : 2 : 3 ratios and allowing us to conjecture a formula for the probabilities for the general N-rows case (which will be discussed and obtained next.)

Now you may be wondering if there's some way to automatically get the exact probability matrix in the limit for a given number N of rows (for reasonable N and running times,) and indeed there is a simple procedure we might try out. First create a copy of my program and edit these five lines to be as follows:

10  DESTROY ALL @ OPTION BASE 1 @ INPUT "Rows=";M @ DIM A(M,M),B(M,M)
15  MAT A=(2/(M*(M+1))) @ W=M-1 @ K=1E-6 @ FOR I=1 TO INF @ MAT B=ZER

70  NEXT Y @ NEXT X @ MAT A=A-B @ DIM A(M*M) @ DISP I;CNORM(A) @ IF RES<K THEN 80
75  DIM A(M,M) @ MAT A=B @ NEXT I
80  FOR I=1 TO M @ FOR J=1 TO I @ DISP FRAC$(B(I,J),6);" "; @ NEXT J @ DISP @ NEXT I

When run, the program asks for the number of rows N and then initializes the matrix probabilities to be the same for all grid locations at the very beginning (and of course adding up to 1) since the limit probability matrix after infinite steps clearly does not depend on the starting position(s) and initially assuming uniformly distributed probabilities greatly speeds up the convergence and accuracy.

Once the initialization is over the program then computes the probability matrix for steps 1, 2, 3, ...., comparing each matrix with the previous one. When the difference is less than a hardcoded tolerance (K=1E-6 at line 15) the process is over and the limit probability matrix is output in exact rational form.

Note: if the rational matrix displayed doesn't look correct (the probabilities don't comply with the 1:2:3 ratios) you can fine-tune the tolerance (say K=1E-7 or smaller, the running time will possibly increase) and/or the accuracy in the conversion to rational form (change the parameter 6 in FRAC$ to some other value, say 7) and run the program again.

While it runs, the program will display the step number and the current difference after each step so you can see it converging to zero, and once it meets the tolerance (will take a long while for large enough N) it will display the limit probability matrix in exact rational form. Let's run it for N=15 rows in FIX 6:

>FIX 6 @ RUN


Rows=15

Step        Difference
----------------------
1.000000     0.991667
2.000000     0.043056
3.000000     0.020833
4.000000     0.012770
...
128.000000   0.000001
129.000000   0.000001

Limit Probability matrix:

1/315
2/315   2/315
2/315   1/105   2/135
    ...
2/315   1/105 ...   1/105   2/315
1/315   2/315 ... ... ...   2/315   1/315

and we see that it took 130 steps (not infinite !) to achieve the specified tolerance. The displayed rational limit probability matrix is nevertheless exact.

Now that we have a working program, we can create a simplified version which just checks the difference of the probability for only the top corner location (1,1) at successive steps, simply comparing A(1,1) vs. B(1,1) instead of considering the whole matrices, and once the tolerance is met we simply output B(1,1), the top corner's probability. If we run this simpler and faster program for various numbers of rows, we get:

N rows    5      7      9       10      11      12      13      14      15
------------------------------------------------------------------------------
P(1,1)    1/30   1/63   1/108   1/135   1/165   1/198   1/234   1/273   1/315

and a little experimentation quickly reveals a pattern for the denominators, e.g.:

P(N=9) -> 108 = 9*24/2,   P(N=10) -> 135 = 10*27/2,   P(N=11) -> 165 = 11*30/2

so we conjecture that the probability in the limit for the top corner of an N-row grid is:

P = 2/(3*N*(N-1))

and taking into account the previously stablished 1:2:3 ratio we finally get

P(corners) = 2/(3*N*(N-1)),   P(edges) = 2*P(corners),   P(inner) = 3*P(corners)

and finally we can create a much simpler and faster program (5 lines or less in all) which will accept N and non-iteratively proceed to immediately display the corresponding limit probability matrix. Checking it for the N=30 rows case we get these probabilities:

P(corners) = 1/1305,   P(edges) = 2/1305  and   P(inner) = 3/1305 = 1/435

which our final program computes and uses to fill up and output the exact probability matrix for the 30-row grid in very little time. Doing the same for a 100,000-row grid would be equally fast.


Well, I hope this has provided a good example of how you can use your vintage HP calc to do some sleuthing and get nice symbolic results in the spirit of Experimental Mathematics. Once you get the numeric results, conjecturing the symbolic ones and afterwards attempting to prove them is that much easier.  

Regards.
V.