----------------------------------------------------------------------
--- Knud van Eeden --- 21 February 2008 - 10:09 pm -------------------
BBCBASIC: Windows: Mathematics: Distribution: Binomial: Hexstat [random walk on a straight line]
---
Given:
The following question was asked by George Dimitriadis:
To simulate the binomial distribution, I used the following code.
It mimics the effects of a HEXSTAT, where an object falling has two
possible choices; left or right. I have used 2*rnd(2)-3 to represent
these choices. However, the shape obtained always appears to be
slightly negatively skewed. I assumed that it will usually appear to be
symmetrical. Is my code invalid?
--- cut here: begin --------------------------------------------------
mode2:vdu5
dim x(15000),y(15000)
for i=0 to 12000
x(i)=0:y(i)= 0
nexti
repeat
x0=800:y0=1000
repeat
r=2*rnd(2)-3
x0=x0+r:y0=y0- .5
until y0<300
y(x0)=y(x0)+ 1
rectangle fill x0-150,50,75, y(x0)
until 0
--- cut here: end ----------------------------------------------------
---
Explanation:
To give a historic overview to support comprehension and the problem
solving process:
This started with a question about a program using the hexstat, if some
calculations were maybe wrong:
It turned out that using RND() led to asymmetric and translated
results, where the result should have been symmetric. So better use
Mersenne Twister or possibly the exact solution.
At first I assumed it very unlikely that something was wrong with
RND(), so assumed instead something else.
E.g. something was missing in the program(s) used.
I analyzed his program (which calculated the distribution of the balls
falling through the hexstat) and realized that calculating where a ball
ends in that hexstat was the same as a man walking on a straight line
and randomly making 1 step backwards or forwards and then noting where
he ends after a given amount of these steps. -It is the same because a
ball falling down in the hexstat can always only go to the left or to
the right from where it is. -And that is the same as the man on the
straight line, who can also always only go to the left or to the right
from where he is.
I calculated a bit and saw zeroes (0) occurring in the end result. That
indicated that there were positions on that straight line where the man
never could arrive. A bit strange, so I thought something has been left
out maybe in my program (maybe something with sum of probabilities not
equal to one or something). I thought that can not be correct, because
you can reach each position on the straight line. Well, actually the
straight line is just part of the story, it only shows the movement in
one direction (but is a good model to get simple calculations, so stick
to it). The hexstat ball is actually moving in 2 directions, to the
left or right and also from top to bottom. So a better model is using a
man hopping on one leg, where he randomly hops to the left or to the
right, from one side of a rectangular field to the other. Now at any
point in a rectangular field, you can hop to the left stone, forward
stone and right stone. As you do not allow to hop to the forward stone,
this explains the zeroes in the result. Or thus expressed in
probabilities you have made the choice (1/2, 0, 1/2). Note that the sum of
the probabilities has to be one, in other words you are taken into
account all possibilities. Another choice would be that you are also
allowed to hop forward. Then the probabilities would be (1/3, 1/3,
1/3), and in that case no zeroes will result, in otherwords you can
reach all positions on the straight line (within the bounds stepmax
to the left to stepmax to the right, plus one for your starting position).
The straight line is the basis of the graphs below. At each step
position on that straight line you indicate with a bar above it which
percentage of the total amount of steps reached that position. You see
that the more steps, the more the graph will get pressed together
towards the center.
Then I searched on the Internet for more about hexstat and this random
walk.
And saw on the Wolfram webpage actually a graph of the random walk.
It turned out to be a Pascal triangle interspersed with zeroes and
divided by 2^(steps).
I calculated it using RND() and got this asymmetric results again and
again. There could not really be something wrong with the calculations
as far as I could tell. The question of George, and the persistence of
Barry Newton, Michael Hutton and David finally at the end convinced me
that it was the RND() function used that was the culprit. There seems
not much there can be done about it (I am checking now with the
original Acorn BBCBASIC micro, but it is already running for 2 days and
nights in a row to calculate 100000 trials, and still nothing coming
up, I will see if some different result turns up there), as the
original BBCBASIC RND() is used, and it has to be kept for backwards
compatibility with legacy BBCBASIC programs).
So better to use the Mersenne Twister random generator. It is slower,
but the version used could in a reasonable time (something like 30
minutes) calculate e.g. 100000 trials, and looks quite OK. Maybe some
time it will be built in in BBCBASIC or written completely in assembler
for speed purposes.
For comparison purposes I created the exact solution (currently maximal
170 steps, but if you rewrite the used functions it could certainly be
increased, e.g. by using the Sterling approximation (tried this, but
turns out to be accepting N=143 maximal against N=170 for the usual
implementation, so skipping this) for the faculty
instead).
For comparison purposes I also created it in other programming
languages like Delphi (very fast, shows the expected graph after pressing
a few times the button) and PHP (and possibly will implement it
in Java, C++, Perl and Ruby also).
---
Calculations:
If you use the expression
x0 + r
you get something like a random walk, on a straight line, one step to
the left and one step to the right, from where you are at that moment.
===
If you checked 'hexstat' structure
E.g.
http://www.mip.berkeley.edu/physics/C+55+05.html
you see indeed that the behavior can be described as a random walk on a straight line,
starting from a fixed point (e.g. x0),
where you each time can take 1 step to the left,
or 1 step to the right.
After N such decisions or thus after N steps, thus e.g. on the Nth row,
splitted out vertically, you have your hexstat.
Or otherwise, if you squeeze the hexstat on a straight line
(or thus project the step vectors on the x-axis, see drawing),
you have your random walk.
Thus
x0 + r
or thus (the next position is the old position minus or plus one step)
where
r=-1
or
r=+1
and r randomly chosen
is a correct representation of the model.
===
To confirm the hypothesis that the zeroes introduced (thus indicating
places where you can not arrive, using the constraints of the steps
used) are correct, I did a search on Google for
"binomial distribution" "random walk"
and see e.g.
http://mathworld.wolfram.com/RandomWalk1-Dimensional.html
for the same graph with zeroes inbetween.
---
steps -5 -4 -3 -2 -1 0 1 2 3 4 5
0 1
1 1/2 0/2 1/2
2 1/4 0/2 2/4 0/2 1/4
3 1/8 0/8 3/8 0/8 3/8 0/8 1/8
4 1/16 0/16 4/16 0/16 6/16 0/16 4/16 0/16 1/16
5 1/32 0/32 5/32 0/32 10/32 0/32 10/32 0/32 5/32 0/32 1/32
===
Here e.g. row 5 is calculated from row 4, working from the left to the right,
by multiplying each value with 1/2, 0, 1/2 respectively, and putting the result
respectively left, below and right. Finally adding the sum vertically
4 1/16 0/16 4/16 0/16 6/16 0/16 4/16 0/16 1/16
| | | | | | | | | | |
| | | | | | | | | | |
v v v | | | | | | | |
+--------------+--+ | | | | | | | |
5 |1/32 0/32 1/32| | | | | | | | |
+--------------+--+ | | | | | | | |
1/32 | | | | | | | | | |
v v | | | | | | | |
+----------------+--+ | | | | | | |
|0/32 0/32 0/32| | | | | | | |
+----------------+--+ | | | | | | |
0/32 | | | | | | | | |
v v | | | | | | |
+---------------+--+ | | | | | |
|4/32 0/32 4/32| | | | | | |
+---------------+--+ | | | | | |
5/32 | | | | | | | |
v v | | | | | |
+---------------+--+ | | | | |
|0/32 0/32 0/32| | | | | |
+---------------+--+ | | | | |
0/32 | | | | | | |
v v | | | | |
+----------------+--+ | | | |
|6/32 0/32 6/32| | | | |
+----------------+--+ | | | |
10/32 | | | | | |
v v | | | |
+---------------+--+ | | |
|0/32 0/32 0/32| | | |
+---------------+--+ | | |
0/32 | | | | |
v v | | |
+--------------+--+ | |
|4/32 0/32 4/32| | |
+--------------+--+ | |
10/32 | | | |
v v | |
+-------------+--+ |
|0/32 0/32 0/32| |
+-------------+--+ |
0/32 | | |
v v v
+-----------------+
|1/32 0/32 1/32|
+-----------------+
5/32 0/32 1/32
===
Note:
In this particular case, thus for probability ( 1/2, 0, 1/2 ), you can also
build the graph by creating the Pascal triangle
1
1 1
1 2 1
1 3 3 1
1 4 6 4 1
1 5 10 10 5 1
Putting zeroes in between
1
1 0 1
1 0 2 0 1
1 0 3 0 3 0 1
1 0 4 0 6 0 4 0 1
1 0 5 0 10 0 10 5 0 1
and dividing each row 1 by 1/2, row 2 by (1/2).(1/2), row 3 by (1/2).(1/2).(1/2),
... row N by (1/2)^N
1
1 0 1
- -
2 2
1 0 2 0 1
- - -
4 4 4
1 0 3 0 3 0 1
- - - -
8 8 8 8
1 0 4 0 6 0 4 0 1
-- -- -- -- --
16 16 16 16 16
1 0 5 0 10 0 10 5 0 1
-- -- -- -- -- --
32 32 32 32 32 32
...
---
To conclude, this is thus a simplest program for the 2 steps variation,
tested for row 5 in the Wolfram graph.
--- cut here: begin --------------------------------------------------
stepMinI = 1
stepMaxI = 5
:
REM this 1-dimensional array represents thus the straight line. You note where you arrive after stepsMaxI steps.
DIM x( stepMaxI + 1 + stepMaxI ) : REM you can not go more than stepMaxI steps to the left (all -1), similar not more than stepsMaxI to the right (all +1)
:
trialMinI = 1
trialMaxI = 100000
:
FOR trialI = trialMinI TO trialMaxI
x0 = stepMaxI + 1 : REM Put x0 in the middle, so keep x0 (used as array index here) positive, so that array index is not negative
FOR stepI = stepMinI TO stepMaxI
r = RND( 2 )
IF ( r = 1 ) THEN x0 += -1 ELSE x0 += +1
NEXT stepI
x( x0 ) += 1 : REM increase counter where you arrive after stepsMaxI steps
NEXT trialI
:
FOR stepI = stepMinI TO ( stepMaxI + 1 + stepMaxI )
PRINT stepI; " "; x( stepI )
NEXT stepI
:
sum = 0
FOR stepI = stepMinI TO ( stepMaxI + 1 + stepMaxI )
sum += x( stepI )
NEXT stepI
PRINT "sum = "; sum
:
FOR stepI = stepMinI TO ( stepMaxI + 1 + stepMaxI )
PRINT x( stepI ) / sum
NEXT stepI
END
--- cut here: end ----------------------------------------------------
===
If you run the program, it gives e.g.
0.03136
0
0.15502
0
0.31109
0
0.31422
0
0.15711
0
0.0312
---
This is the program which is analogue to the zero graph at Wolfram.
Run e.g. the program, and compare the result for row 5.
That gives 1/32, 0, 5/32, 0, 10/32, 0, 10/32, 0, 5/32/, 0, 1/32
PRINT 1/32 gives 0.03125
PRINT 5/32 gives 0.15625
PRINT 10/32 gives 0.3125
Thus what you expect theoretically is in the long run is
0.03125
0
0.15625
0
0.3125
0
0.3125
0
0.15625
0
0.03125
===
Another program drawing a graph (see also drawing below).
Note: this graph is supposed to be symmetric (it is thus a Pascal triangle,
interspersed with zeroes and scaled), and the red line should be in the middle.
When tested in BBCBASIC v2.0 and BBCBASIC v5.80a, in both cases the
graph is translated away from the middle, and is asymmetric.
This behavior is due to the RND() random generator used in
BBCBASIC.
When using this built-in RND() its results should therefore be
further analyzed and interpreted very carefully.
===
--- cut here: begin --------------------------------------------------
MODE8
:
trialMinI = 1
trialMaxI = 1000000
:
stepMinI = 1
stepMaxI = 640
:
yMaxI = -999999
yScreenMaxI = 800
:
DIM x( stepMaxI + 1 + stepMaxI )
:
FOR trialI = trialMinI TO trialMaxI
x0 = stepMaxI + 1
FOR stepI = stepMinI TO stepMaxI
CASE RND( 2 ) OF
WHEN 1
x0 += 1
WHEN 2
x0 -= 1
ENDCASE
NEXT
x( x0 ) += 1
NEXT
:
FOR stepI = stepMinI TO ( stepMaxI + 1 + stepMaxI )
IF yMaxI < x( stepI ) THEN yMaxI = x( stepI )
NEXT
:
FOR stepI = stepMinI TO ( stepMaxI + 1 + stepMaxI )
MOVE stepI, 0
DRAW stepI, ( x( stepI ) / yMaxI ) * yScreenMaxI
NEXT
:
GCOL1
MOVE stepMaxI + 1, 0
DRAW stepMaxI + 1, 800
:
--- cut here: end ----------------------------------------------------
===
Another program drawing a graph (see also drawing below).
The same restrictions occur in the graph of this program (because of
use of built-in RND()).
===
--- cut here: begin --------------------------------------------------
MODE8
:
trialMinI = 1
trialMaxI = 100000
:
stepMinI = 1
stepMaxI = 640
:
yMaxI = -999999
yScreenMaxI = 800
:
DIM x( stepMaxI + 1 + stepMaxI )
:
FOR trialI = trialMinI TO trialMaxI
x0 = stepMaxI + 1
FOR stepI = stepMinI TO stepMaxI
IF RND( 1 ) > 0.5 THEN
x0 += 1
ELSE
x0 -= 1
ENDIF
NEXT stepI
x( x0 ) += 1
NEXT
:
FOR stepI = stepMinI TO ( stepMaxI + 1 + stepMaxI )
IF yMaxI < x( stepI ) THEN yMaxI = x( stepI )
NEXT
:
FOR stepI = stepMinI TO ( stepMaxI + 1 + stepMaxI )
MOVE stepI, 0
DRAW stepI, ( x( stepI ) / yMaxI ) * yScreenMaxI
NEXT
:
GCOL1
MOVE stepMaxI + 1, 0
DRAW stepMaxI + 1, 800
:
--- cut here: end ----------------------------------------------------
===
Another program drawing a graph (see also drawing below).
This program uses a Mersenne Twister random generator
(this routine is much more accurate, but slower,
as not written completely in assembler,
as the built-in RND() function is)
===
--- cut here: begin --------------------------------------------------
MODE8
:
seedI = FNMathGetRandomBetween1AndNMersenneTwisterR( -ABS(TIME) )
:
trialMinI = 1
trialMaxI = 10000
:
stepMinI = 1
stepMaxI = 640
:
yMaxI = -999999
yScreenMaxI = 800
:
DIM x( stepMaxI + 1 + stepMaxI )
:
FOR trialI = trialMinI TO trialMaxI
x0 = stepMaxI + 1
FOR stepI = stepMinI TO stepMaxI
IF ( FNMathGetRandomBetween1AndNMersenneTwisterR( 2 ) == 1 ) THEN
x0 += 1
ELSE
x0 -= 1
ENDIF
NEXT stepI
x( x0 ) += 1
NEXT
:
FOR stepI = stepMinI TO ( stepMaxI + 1 + stepMaxI )
IF yMaxI < x( stepI ) THEN yMaxI = x( stepI )
NEXT
:
FOR stepI = stepMinI TO ( stepMaxI + 1 + stepMaxI )
MOVE stepI, 0
DRAW stepI, ( x( stepI ) / yMaxI ) * yScreenMaxI
NEXT
:
GCOL1
MOVE stepMaxI + 1, 0
DRAW stepMaxI + 1, 800
:
END
:
:
:
REM library: math: get: random: between1: and: n: mersenne: twister (filenamemacro=getmatt.bbc) <author>Richard Russell</author> [kn, ri, mo, 25-02-2008 02:24:52]
DEF FNMathGetRandomBetween1AndNMersenneTwisterR( R% )
REM e.g. seed% = FNMathGetRandomBetween1AndNMersenneTwisterR( -ABS(TIME) )
REM e.g. FOR N% = 1 TO 30
REM e.g. PRINT FNMathGetRandomBetween1AndNMersenneTwisterR( 2 )
REM e.g. NEXT
REM e.g. END
PRIVATE MT%(), J%
LOCAL I%, Y%
DIM MT%(623)
IF R%<0 THEN
LOCAL A%, B%, C%, P%, Q%
DIM P% LOCAL 20
[OPT 2:.Q% mov edx,eax:shr eax,30:xor eax,edx:mul ebx:add eax,ecx:ret:]
MT%(0) = -R%
B% = 1812433253
FOR C% = 1 TO 623
A% = MT%(C%-1)
MT%(C%) = USR(Q%)
NEXT
J% = 624
= R%
ENDIF
REPEAT
IF J% = 624 THEN
REM Generate an array of 624 untempered numbers:
FOR I% = 0 TO 623
Y% = (MT%(I%) AND &80000000) OR (MT%((I%+1) MOD 624) AND &7FFFFFFF)
MT%(I%) = MT%((I% + 397) MOD 624) EOR (Y% >>> 1)
IF (Y% AND 1) MT%(I%) EOR= &9908B0DF
NEXT
J% = 0
ENDIF
REM Extract a tempered pseudorandom 32-bit number:
Y% = MT%(J%)
J% += 1
Y% EOR= Y% >>> 11
Y% EOR= (Y% << 7) AND &9D2C5680
Y% EOR= (Y% << 15) AND &EFC60000
Y% EOR= Y% >>> 18
IF R%<=0 THEN = Y% : REM Return 32-bit integer
UNTIL (Y% AND &7FFFFFFF) < ((&7FFFFFFF DIV R%) * R%)
= (Y% AND &7FFFFFFF) MOD R% + 1 : REM Return value in range 1 to R%
:
--- cut here: end ----------------------------------------------------
---
You see thus that the '0' is like the pin in a pinball machine (which
shows clearly in the Wolfram graph), where it bounces against, and can
go in 2 directions, left or right.
Or in general, the '0' are the pins.
0
0 0
0 0 0
0 0 0 0
To see this better, imagine that you build your hexstat using
nails.
Throwing e.g. 'x' through it will give e.g. the result
0
0 0
0 0 0
0 x 0 x 0 x 0
You see that the nail is thus the '0', that is a position where the 'x'
can not come.
===
To illustrate the choice RND(2), or thus -1 and +1 steps.
Similar to hop on a rectangular field (e.g. young people jumping on
one leg), when you are not allowed to hop straigt forward, but only right forward.
or left forward.
You can only hop in 2 directions.
Thus a probability of 1/2 to go left, a probability of 0 to go forward, a
probability of 1/2 to go right (note that 1/2 + 0 + 1/2 = 1, thus a
total probability of 1).
You see thus clearly there are certain positions which you can not reach
(usually the straight forward hop).
===
0 hops
+-----+-----+-----+-----+-----+-----+-----+
| | | | | | | |
| | | |start| | | |
| | | | | | | |
+-----+-----+-----+-----+-----+-----+-----+
| | | | | | | |
| | | | | | | |
| | | | | | | |
+-----+-----+-----+-----+-----+-----+-----+
| | | | | | | |
| | | | | | | |
| | | | | | | |
+-----+-----+-----+-----+-----+-----+-----+
| | | | | | | |
| | | | | | | |
| | | | | | | |
+-----+-----+-----+-----+-----+-----+-----+
1 hop
+-----+-----+-----+-----+-----+-----+-----+
| | | | | | | |
| | | |start| | | |
| | | | | | | |
+-----+-----+-----+-----+-----+-----+-----+
| | | | | | | |
| | |left | |right| | |
| | | | | | | |
+-----+-----+-----+-----+-----+-----+-----+
| | | | | | | |
| | | | | | | |
| | | | | | | |
+-----+-----+-----+-----+-----+-----+-----+
| | | | | | | |
| | | | | | | |
| | | | | | | |
| | | | | | | |
+-----+-----+-----+-----+-----+-----+-----+
2 hops
+-----+-----+-----+-----+-----+-----+-----+
| | | | | | | |
| | | |start| | | |
| | | | | | | |
+-----+-----+-----+-----+-----+-----+-----+
| | | | | | | |
| | |left | |right| | |
| | | | | | | |
+-----+-----+-----+-----+-----+-----+-----+
| | | |left | | | |
| |left | |right| |right| |
| | | | | | | |
+-----+-----+-----+-----+-----+-----+-----+
| | | | | | | |
| | | | | | | |
| | | | | | | |
+-----+-----+-----+-----+-----+-----+-----+
3 hops
+-----+-----+-----+-----+-----+-----+-----+
| | | | | | | |
| | | |start| | | |
| | | | | | | |
+-----+-----+-----+-----+-----+-----+-----+
| | | | | | | |
| | |left | |right| | |
| | | | | | | |
+-----+-----+-----+-----+-----+-----+-----+
| | | |left | | | |
| |left | |right| |right| |
| | | | | | | |
+-----+-----+-----+-----+-----+-----+-----+
| | |left | |right| | |
| | |left | |left | | |
|left | |right| |right| |right|
| | |right| |left | | |
+-----+-----+-----+-----+-----+-----+-----+
and so on...
===
Exact calculation:
Pascal triangle
1
1 1
1 2 1
1 3 3 1
1 4 6 4 1
1 5 10 10 5 1
1 6 15 20 15 6 1
...
===
Pascal triangle using combinatorics C( m, n )
There are several methods:
===
METHOD 1: CALCULATE THE INDIVIDUAL COEFFICIENTS OF THE PASCAL TRIANGLE
using
C( n, m ) = n! / m! / ( n - m )!
you can write the values in the Pascal triangle then as follows:
C(0,0)
C(1,0) C(1,1)
C(2,0) C(2,1) C(2,2)
C(3,0) C(3,1) C(3,2) C(3,3)
C(4,0) C(4,1) C(4,2) C(4,3) C(4,4)
C(5,0) C(5,1) C(5,2) C(5,3) C(5,4) C(5,5)
C(6,0) C(6,1) C(6,2) C(6,3) C(6,4) C(6,5) C(6,6)
...
C(N,0) C(N,1) C(N,2) C(N,3) ... C(N,N-2) C(N,N-1) C(N,N)
===
Using a for next loop to show the values in a particular horizontal row:
FOR I% = 0 TO N%
PRINT FNMathGetCombinationI( N%, I% )
NEXT I%
===
Pascal triangle using combinatorics C( m, n ) and interspersed with zeroes
C(0,0)
C(1,0) 0 C(1,1)
C(2,0) 0 C(2,1) 0 C(2,2)
C(3,0) 0 C(3,1) 0 C(3,2) 0 C(3,3)
C(4,0) 0 C(4,1) 0 C(4,2) 0 C(4,3) 0 C(4,4)
C(5,0) 0 C(5,1) 0 C(5,2) 0 C(5,3) 0 C(5,4) 0 C(5,5)
C(6,0) 0 C(6,1) 0 C(6,2) 0 C(6,3) 0 C(6,4) 0 C(6,5) 0 C(6,6)
...
C(N,0) 0 C(N,1) 0 C(N,2) 0 C(N,3) 0 ... 0 C(N,N-2) 0 C(N,N-1) 0 C(N,N)
===
Using a for next loop to show the values in a particular horizontal row:
FOR I% = 0 TO N%
:
PRINT FNMathGetCombinationI( N%, I% )
:
REM if not the last term, the next value is a zero
IF ( I% <> N% ) THEN
PRINT 0
ENDIF
:
NEXT I%
===
Pascal triangle using combinatorics C( m, n ) and interspersed with zeroes and divided by power of 2
C(0,0) / 2^0 : step 0
C(1,0) 0 C(1,1) / 2^1 : step 1
C(2,0) 0 C(2,1) 0 C(2,2) / 2^2 : step 2
C(3,0) 0 C(3,1) 0 C(3,2) 0 C(3,3) / 2^3 : step 3
C(4,0) 0 C(4,1) 0 C(4,2) 0 C(4,3) 0 C(4,4) / 2^4 : step 4
C(5,0) 0 C(5,1) 0 C(5,2) 0 C(5,3) 0 C(5,4) 0 C(5,5) / 2^5 : step 5
C(6,0) 0 C(6,1) 0 C(6,2) 0 C(6,3) 0 C(6,4) 0 C(6,5) 0 C(6,6) / 2^6 : step 6
... ...
C(N,0) 0 C(N,1) 0 C(N,2) 0 C(N,3) 0 ... 0 C(N,N-2) 0 C(N,N-1) 0 C(N,N) / 2^N : step N
===
Using a for next loop to show the values in a particular horizontal row:
--- cut here: begin --------------------------------------------------
FOR I% = 0 TO N%
:
PRINT FNMathGetCombinationI( N%, I% ) / 2^N%
:
REM if not the last term, the next value is a zero
IF ( I% <> N% ) THEN
PRINT 0
ENDIF
:
NEXT I%
--- cut here: end ----------------------------------------------------
===
E.g. if 640 steps then N% = 640
and you get
--- cut here: begin --------------------------------------------------
FOR I% = 0 TO 640
:
PRINT FNMathGetCombinationI( 640, I% ) / 2^640
:
REM if not the last term, the next value is a zero
IF ( I% <> 640 ) THEN
PRINT 0
ENDIF
:
NEXT I%
--- cut here: end ----------------------------------------------------
Because of this large numbers (like 2^N%), this will usually limit the
maximum amount of steps which you can do (unless you e.g. use or create
special data structures or use symbolic mathematics programs).
E.g. in a symbolic mathematics program like Maple you can write
--- cut here: begin --------------------------------------------------
for i from 0 to 640 do print( ( C( 640, i ) / 2^640 ) * 1.0 ); print ( 0 ) ; od;
--- cut here: end ----------------------------------------------------
This gives at the extremities very small values near zero, with at the
center the following values (use this e.g. for checking your results).
--- cut here: begin --------------------------------------------------
...
-7
.8145171533 10
0
-6
.1217021739 10
0
-6
.1806664054 10
0
-6
.2664655091 10
0
-6
.3904744577 10
0
-6
.5685068733 10
0
-6
.8223820800 10
0
-5
.1181978807 10
0
-5
.1687901554 10
0
-5
.2394909373 10
0
-5
.3376282012 10
0
-5
.4729323868 10
0
-5
.6582230608 10
0
-5
.9102564259 10
0
.00001250759756
0
.00001707679371
0
.00002316667971
0
.00003122834481
0
.00004182774651
0
.00005566892808
0
.00007362014039
0
.00009674271156
0
.0001263223176
0
.0001639020752
0
.0002113166041
0
.0002707258984
0
.0003446475090
0
.0004359851881
0
.0005480518033
0
.0006845840071
0
.0008497458829
0
.001048118615
0
.001284673164
0
.001564723023
0
.001893854417
0
.002277831773
0
.002722477016
0
.003233522191
0
.003816436055
0
.004476226696
0
.005217223682
0
.006042844938
0
.006955355078
0
.007955623536
0
.009042892085
0
.01021456249
0
.01146601551
0
.01279047274
0
.01417891222
0
.01562004756
0
.01710037887
0
.01860432098
0
.02011441197
0
.02161160121
0
.02307561291
0
.02448537703
0
.02581951617
0
.02705687318
0
.02817706219
0
.02916102309
0
.02999155857
0
.03065383272
0
.03113581123
0
.03142862451
0
.03152683896
0
.03142862451
0
.03113581123 <- here is the maximum
0
.03065383272
0
.02999155857
0
.02916102309
0
.02817706219
0
.02705687318
0
.02581951617
0
.02448537703
0
.02307561291
0
.02161160121
0
.02011441197
0
.01860432098
0
.01710037887
0
.01562004756
0
.01417891222
0
.01279047274
0
.01146601551
0
.01021456249
0
.009042892085
0
.007955623536
0
.006955355078
0
.006042844938
0
.005217223682
0
.004476226696
0
.003816436055
0
.003233522191
0
.002722477016
0
.002277831773
0
.001893854417
0
.001564723023
0
.001284673164
0
.001048118615
0
.0008497458829
0
.0006845840071
0
.0005480518033
0
.0004359851881
0
.0003446475090
0
.0002707258984
0
.0002113166041
0
.0001639020752
0
.0001263223176
0
.00009674271156
0
.00007362014039
0
.00005566892808
0
.00004182774651
0
.00003122834481
0
.00002316667971
0
.00001707679371
0
.00001250759756
0
-5
.9102564259 10
0
-5
.6582230608 10
0
-5
.4729323868 10
0
-5
.3376282012 10
0
-5
.2394909373 10
0
-5
.1687901554 10
0
-5
.1181978807 10
0
-6
.8223820800 10
0
-6
.5685068733 10
0
-6
.3904744577 10
0
-6
.2664655091 10
0
-6
.1806664054 10
0
-6
.1217021739 10
0
-7
.8145171533 10
...
--- cut here: end ----------------------------------------------------
===
A simplest program which calculates the exact values is thus:
--- cut here: begin --------------------------------------------------
INPUT "steps total (maximum to choose is 34) = "; stepMaxI
:
powerOfTwoI = 2^stepMaxI
:
FOR stepI = 0 TO stepMaxI
:
PRINT; FNMathGetCombinationI( stepMaxI, stepI ) / powerOfTwoI
:
REM if not the last term, the next value is a zero
IF ( stepI <> stepMaxI ) THEN
PRINT; 0
ENDIF
:
NEXT stepI
:
END
:
:
:
DEF FNMathGetCombinationI( N%, P% )
= FNMathGetFacultyI( N% ) / FNMathGetFacultyI( P% ) / FNMathGetFacultyI( N% - P% )
:
DEF FNMathGetFacultyI( I% )
IF ( I% < 0 ) THEN PRINT "please enter only 0 or positive integers" : = 0
IF ( I% = 0 ) THEN = 1
IF ( I% = 1 ) THEN = 1
= I% * FNMathGetFacultyI( I% - 1 )
:
--- cut here: end ----------------------------------------------------
===
A simple program which calculates the exact values and shows the result
in a graph (by using double precision (*FLOAT 64) the total amount of
steps is increased)
Another program drawing a graph (see also drawing below).
--- cut here: begin --------------------------------------------------
*FLOAT 64
:
MODE 8
:
INPUT "steps total (maximum to choose is 170) = "; stepMaxI
:
powerOfTwoI = 2^stepMaxI
:
xScreenMinI = 0
xScreenMaxI = 1000
yScreenMinI = 0
yScreenMaxI = 800
:
stepMinI = 1
:
I = stepMinI - 1 : REM let the array counter start at stepMinI
:
yMaxI = -999999
:
DIM x( stepMaxI + 1 + stepMaxI )
:
FOR stepI = 0 TO stepMaxI
:
I = I + 1
x( I ) = FNMathGetCombinationI( stepMaxI, stepI ) / powerOfTwoI
:
REM if not the last term, the next value is a zero
IF ( stepI <> stepMaxI ) THEN
I = I + 1
x( I ) = 0
ENDIF
:
NEXT stepI
:
FOR stepI = stepMinI TO stepMaxI + 1 + stepMaxI
IF yMaxI < x( stepI ) THEN
yMaxI = x( stepI )
ENDIF
NEXT stepI
:
FOR stepI = stepMinI TO ( stepMaxI + 1 + stepMaxI )
MOVE FNx( stepI ), FNy( yScreenMinI )
DRAW FNx( stepI ), FNy( x( stepI ) )
NEXT stepI
:
GCOL 1
MOVE FNx( stepMaxI + 1 ), 0
DRAW FNx( stepMaxI + 1 ), yScreenMaxI
:
END
:
:
:
DEF FNx( x )
= x / ( stepMaxI + 1 + stepMaxI ) * xScreenMaxI
:
DEF FNy( y )
= ( y / yMaxI ) * yScreenMaxI
:
DEF FNMathGetCombinationI( N%, P% )
= FNMathGetFacultyI( N% ) / FNMathGetFacultyI( P% ) / FNMathGetFacultyI( N% - P% )
:
DEF FNMathGetFacultyI( I% )
IF ( I% < 0 ) THEN PRINT "please enter only 0 or positive integers" : = 0
IF ( I% = 0 ) THEN = 1
IF ( I% = 1 ) THEN = 1
= I% * FNMathGetFacultyI( I% - 1 )
:
--- cut here: end ----------------------------------------------------
===
METHOD 2: BUILD THE WHOLE PASCAL TRIANGLE IN AN ARRAY
You can build the Pascal triangle in an array, by starting with 1, 1
and 1, and then continuously adding, as usual.
The idea is to add adjacent numbers, and put that sum in the row, below
the last of that two numbers.
Ideas:
1. You input the total number of rows n
2. Use a matrix of n rows and n columns
(this matrix should be filled with only zeros when you start)
3. In the first column all 1s.
4. Use the following little program (the idea is to add adjacent numbers, and put that sum in the row, below the last of that two numbers)
--- cut here: begin --------------------------------------------------
rowMaxI = 4
n = rowMaxI + 1
DIM matrix( n, n )
FOR row = 1 TO n
matrix( row, 1 ) = 1
FOR column = 2 TO row
matrix( row, column ) = matrix( row - 1, column - 1 ) + matrix( row - 1, column )
NEXT column
NEXT row
:
REM print the result
:
FOR row = 1 TO n
FOR column = 1 TO n
PRINT matrix( row, column ); " ";
NEXT column
PRINT
NEXT row
--- cut here: end ----------------------------------------------------
===
If you run this program you get:
--- cut here: begin --------------------------------------------------
>RUN
1 0 0 0 0
1 1 0 0 0
1 2 1 0 0
1 3 3 1 0
1 4 6 4 1
--- cut here: end ----------------------------------------------------
Using a brute force approach by dimensioning an NxN array
(you store thus the whole Pascal triangle in memory. Actually for the
purposes you only need one or two rows, so a 1xN or 2xN array should
do) you can calculate maximal around 358 rows of the Pascal triangle, or thus
358 steps or thus 358 rows in the hexstat).
It calculates the last row of the stored Pascal triangle, multiplies it
with a power of (1/2), and put zeroes in between.
Another program drawing a graph (see also drawing below).
--- cut here: begin --------------------------------------------------
*FLOAT 64
:
MODE 8
:
INPUT "steps total (maximum to choose is 358) = "; stepMaxI
:
n = stepMaxI + 1
:
xScreenMinI = 0
xScreenMaxI = 1000
yScreenMinI = 0
yScreenMaxI = 800
:
stepMinI = 1
:
I = stepMinI - 1 : REM let the array counter start at stepMinI
:
yMaxI = -999999
:
DIM x( stepMaxI + 1 + stepMaxI )
:
DIM matrix( n, n )
:
FOR rowI = 1 TO n
matrix( rowI, 1 ) = 1
FOR columnI = 2 TO rowI
matrix( rowI, columnI ) = matrix( rowI - 1, columnI - 1 ) + matrix( rowI - 1, columnI )
NEXT columnI
NEXT rowI
:
FOR rowI = n TO n
FOR columnI = 1 TO n
REM Multiplying it with the power of 2 (or thus multiplying it with 1/2, the chosen probability)
:
I = I + 1
x( I ) = matrix( rowI, columnI ) * (1/2)^stepMaxI
:
REM if not the last term, the next value is a zero
IF ( columnI <> n ) THEN
I = I + 1
x( I ) = 0
ENDIF
NEXT columnI
NEXT rowI
:
FOR stepI = stepMinI TO stepMaxI + 1 + stepMaxI
IF yMaxI < x( stepI ) THEN
yMaxI = x( stepI )
ENDIF
NEXT stepI
:
FOR stepI = stepMinI TO ( stepMaxI + 1 + stepMaxI )
MOVE FNx( stepI ), FNy( yScreenMinI )
DRAW FNx( stepI ), FNy( x( stepI ) )
NEXT stepI
:
GCOL 1
MOVE FNx( stepMaxI + 1 ), 0
DRAW FNx( stepMaxI + 1 ), yScreenMaxI
:
END
:
:
:
DEF FNx( x )
= x / ( stepMaxI + 1 + stepMaxI ) * xScreenMaxI
:
DEF FNy( y )
= ( y / yMaxI ) * yScreenMaxI
:
--- cut here: end ----------------------------------------------------
===
METHOD: CALCULATING ONLY THE LAST ROW OF THE PASCAL TRIANGLE
The goal is to increase the total amount of steps you can
calculate. Until now this are 170 and 358 steps maximal.
The following program will be more efficient, by only having
to store the but last row and the last row of the Pascal triangle.
This increases the maximum amount of steps to 1029.
Create an array with 2 rows and N columns.
You start with 1 in the first column in the first row.
Then you use the rule to add the previous row values to the next row.
When finished, you copy the second row back to the first row,
and start again.
Using *FLOAT 64, you can so calculate 1029 rows of the Pascal triangle.
--- cut here: begin --------------------------------------------------
*FLOAT 64
rowMaxI = 1029
:
n = rowMaxI + 1
:
DIM matrix( 2, n )
:
matrix( 1, 1 ) = 1
:
FOR row = 1 TO n
FOR column = 1 TO row
matrix( 2, column ) = matrix( 1, column - 1 ) + matrix( 1, column )
NEXT column
:
FOR column = 2 TO n
matrix( 1, column ) = matrix( 2, column )
NEXT column
NEXT row
:
REM print the result
:
FOR column = 1 TO n
PRINT matrix( 2, column ); " ";
NEXT column
--- cut here: end ----------------------------------------------------
Using this building block, you can create the following program.
Another program drawing a graph (see also drawing below).
--- cut here: begin --------------------------------------------------
*FLOAT 64
:
MODE 8
:
INPUT "steps total (maximum to choose is 1022) = "; stepMaxI
:
n = stepMaxI + 1
:
xScreenMinI = 0
xScreenMaxI = 1000
yScreenMinI = 0
yScreenMaxI = 800
:
stepMinI = 1
:
I = stepMinI - 1 : REM let the array counter start at stepMinI
:
yMaxI = -999999
:
DIM x( stepMaxI + 1 + stepMaxI )
:
DIM matrix( 2, n )
:
matrix( 1, 1 ) = 1
:
FOR rowI = 1 TO n
FOR columnI = 1 TO rowI
matrix( 2, columnI ) = matrix( 1, columnI - 1 ) + matrix( 1, columnI )
NEXT columnI
:
FOR columnI = 2 TO n
matrix( 1, columnI ) = matrix( 2, columnI )
NEXT columnI
NEXT rowI
:
FOR columnI = 1 TO n
:
I = I + 1
REM Multiplying it with the power of 2 (or thus multiplying it with 1/2, the chosen probability)
x( I ) = matrix( 2, columnI ) * (1/2)^stepMaxI
:
REM if not the last term, the next value is a zero
IF ( columnI <> n ) THEN
I = I + 1
x( I ) = 0
ENDIF
:
NEXT columnI
:
FOR stepI = stepMinI TO stepMaxI + 1 + stepMaxI
IF yMaxI < x( stepI ) THEN
yMaxI = x( stepI )
ENDIF
NEXT stepI
:
FOR stepI = stepMinI TO ( stepMaxI + 1 + stepMaxI )
MOVE FNx( stepI ), FNy( yScreenMinI )
DRAW FNx( stepI ), FNy( x( stepI ) )
NEXT stepI
:
GCOL 1
MOVE FNx( stepMaxI + 1 ), 0
DRAW FNx( stepMaxI + 1 ), yScreenMaxI
:
END
:
:
:
DEF FNx( x )
= x / ( stepMaxI + 1 + stepMaxI ) * xScreenMaxI
:
DEF FNy( y )
= ( y / yMaxI ) * yScreenMaxI
:
--- cut here: end ----------------------------------------------------
USING OTHER PROBABILITIES
===
To illustrate the choice RND(3), or thus -1, 0, +1 steps
Similar to hop on a rectangular field (e.g. young people jumping on
one leg), when you are allowed to hop straigt forward, right forward,
or left forward.
You are thus allowed to hop in 3 directions here.
Thus a probability of 1/3 to go left, a probability of 1/3 to go forward, a
probability of 1/3 to go right (note that 1/3 + 1/3 + 1/3 = 1, thus a
total probability of 1).
===
0 hops
+-----+-----+-----+-----+-----+-----+-----+
| | | | | | | |
| | | |start| | | |
| | | | | | | |
+-----+-----+-----+-----+-----+-----+-----+
| | | | | | | |
| | | | | | | |
| | | | | | | |
+-----+-----+-----+-----+-----+-----+-----+
| | | | | | | |
| | | | | | | |
| | | | | | | |
+-----+-----+-----+-----+-----+-----+-----+
| | | | | | | |
| | | | | | | |
| | | | | | | |
+-----+-----+-----+-----+-----+-----+-----+
1 hop
+-----+-----+-----+-----+-----+-----+-----+
| | | | | | | |
| | | |start| | | |
| | | | | | | |
+-----+-----+-----+-----+-----+-----+-----+
| | | | | | | |
| | |left |mid |right| | |
| | | | | | | |
+-----+-----+-----+-----+-----+-----+-----+
| | | | | | | |
| | | | | | | |
| | | | | | | |
+-----+-----+-----+-----+-----+-----+-----+
| | | | | | | |
| | | | | | | |
| | | | | | | |
+-----+-----+-----+-----+-----+-----+-----+
2 hops
+-----+-----+-----+-----+-----+-----+-----+
| | | | | | | |
| | | |start| | | |
| | | | | | | |
+-----+-----+-----+-----+-----+-----+-----+
| | | | | | | |
| | |left |mid |right| | |
| | | | | | | |
+-----+-----+-----+-----+-----+-----+-----+
| | |left |mid |right| | |
| |left |mid |right|mid |right| |
| | | |left | | | |
+-----+-----+-----+-----+-----+-----+-----+
| | | | | | | |
| | | | | | | |
| | | | | | | |
+-----+-----+-----+-----+-----+-----+-----+
and so on...
===
If you use 3 possibilities for the steps: -1, 0 and +1
---
By using the rule that all the time at each point one third goes to the left,
one third goes straight on, one third goes to the right, repeating this for
all the values on the row, and adding the results, working from top to bottom
you can similarly create the Wolfram graph for RND(3).
---
steps -5 -4 -3 -2 -1 0 1 2 3 4 5
0 1
1 1/3 1/3 1/3
2 1/9 2/9 3/9 2/9 1/9
3 1/27 3/27 6/27 7/27 6/27 3/27 1/27
4 1/81 4/81 10/81 16/81 19/81 16/81 10/81 4/81 1/81
5 1/243 5/243 15/243 30/243 45/243 51/243 45/243 30/243 15/243 5/243 1/243
===
To see how to build this graph yourself manually:
Here you see how I calculated e.g. row 5. I all the time multiplied the
value of row 4 with 1/3. And did this 3 times. And wrote the result at
the left, middle and right respectively. I did this for all the
values of row 4. Finally I added the results vertically.
5 1/243 1/243 1/243
4/243 4/243 4/243
10/243 10/243 10/243
16/243 16/243 16/243
19/243 19/243 19/243
16/243 16/243 16/243
10/243 10/243 10/243
4/243 4/243 4/243
1/243 1/243 1/243
---
Here a simplest program:
To include the middle pin, you can choose steps -1, 0, +1, and use
RND(3), and use
r = RND(3)
IF ( r = 1 ) THEN
x0 += -1
ELSE IF ( r = 2 ) THEN
x0 += 0
ELSE
x0 += 1
ENDIF
ENDIF
--- cut here: begin --------------------------------------------------
stepMinI = 1
stepMaxI = 5
:
REM this 1-dimensional array represents thus the straight line. You note where you arrive after stepsMaxI steps.
DIM x( stepMaxI + 1 + stepMaxI ) : REM you can not go more than stepMaxI steps to the left (all -1), similar not more than stepsMaxI to the right (all +1)
:
trialMinI = 1
trialMaxI = 100000
:
FOR trialI = trialMinI TO trialMaxI
x0 = stepMaxI + 1 : REM Put x0 in the middle, so keep x0 (used as array index here) positive, so that array index is not negative
FOR stepI = stepMinI TO stepMaxI
r = RND(3)
IF ( r = 1 ) THEN
x0 += -1
ELSE IF ( r = 2 ) THEN
x0 += 0
ELSE
x0 += 1
ENDIF
ENDIF
NEXT stepI
x( x0 ) += 1 : REM increase counter where you arrive after stepsMaxI steps
NEXT trialI
:
FOR stepI = stepMinI TO ( stepMaxI + 1 + stepMaxI )
PRINT stepI; " "; x( stepI )
NEXT stepI
:
sum = 0
FOR stepI = stepMinI TO ( stepMaxI + 1 + stepMaxI )
sum += x( stepI )
NEXT stepI
PRINT "sum = "; sum
:
FOR stepI = stepMinI TO ( stepMaxI + 1 + stepMaxI )
PRINT x( stepI ) / sum
NEXT stepI
END
--- cut here: end ----------------------------------------------------
===
If you run the program it gives e.g. the values
0.00446
0.02043
0.06137
0.12317
0.18714
0.20839
0.18374
0.12565
0.06061
0.02073
0.00431
compare this e.g. with the theoretical expected values (on row 5),
see row 5 in the Wolfram Graph above.
PRINT 1/243
0.00411522634
PRINT 5/243
0.0205761317
PRINT 15/243
0.0617283951
PRINT 30/243
0.12345679
PRINT 45/243
0.185185185
PRINT 51/243
0.209876543
or thus theoretically if 5 steps, and 3 possibilities for the steps -1, 0 and +1
0.00411522634
0.0205761317
0.0617283951
0.12345679
0.185185185
0.209876543
0.185185185
0.12345679
0.0617283951
0.0205761317
0.00411522634
===
Screencast: see also:
===
Image: see also:
===
Program: see also:
01.exe - RND(2)
02.exe - RND(1)
03.exe - Mersenne Twister
04.exe - Exact: calculated with C(m,n). Maximal 170 steps. (Pascal triangle interspersed with zeroes and divided by power of 2)
05.exe - Exact: calculated by storing whole Pascal triangle in memory and taking last row. Maximal 358 steps. (Pascal triangle interspersed with zeroes and divided by power of 2)
06.exe - Exact: calculating only last row of Pascal triangle. Maximal 1022 steps. (Pascal triangle interspersed with zeroes and divided by power of 2)
===
Video: see also:
===
Internet: see also:
Mathematics: Probability: Distribution: Binomial: Hexstat: Links: Can you give an overview of links?
http://www.knudvaneeden.com/tinyurl.php?urlKey=url000123
----------------------------------------------------------------------