'============================================================================
'   G.H. George                                                 1993 FEB 21
'                                  CLT.BAS            Modified  2000 AUG 10
'   An illustration for the central limit theorem, (many rolls of a loaded
'    die).
'============================================================================

'   Subprograms / Functions:
'          Assign:  set the probabilities for the outcomes of the dice.
'       Histogram:  draw an histogram for the current sample mean  X_n
DECLARE SUB Assign (Mean!)
DECLARE SUB Histogram (n AS INTEGER, Mean!)
DECLARE SUB IntroScreen ()
DECLARE FUNCTION f! (x!, n AS INTEGER)

'   Global variable:
'         P():  array holding probabilities for  n  loaded dice (n = 2^k).
COMMON SHARED P()
DIM P(1 TO 384, 1 TO 7)

'   Global constants:
'        True:  logical true ( = -1)
'       False:  logical false ( = 0)
CONST False = 0, True = NOT False

'   Local variable:
'  SampleSize:  current sample size (doubles each time around the loop).
DIM SampleSize AS INTEGER

CALL Assign(Mean)

LET SampleSize = 1

DO
    CALL Histogram(SampleSize, Mean)
LOOP UNTIL SampleSize > 64

END

SUB Assign (Mean)
'============================================================================
'   Print a welcome message, then:
'   Ask the user to assign probabilities to each of the 6 outcomes of 1 die,
'   then calculate probabilities for each of the  (5n + 1)  outcomes of
'   rolling the die  n  times,  for  n = 1, 2, 4, 8, 16, 32 and 64.
'   Also returns the expected value of  X  for  n = 1 :  Mean = E[X].
'============================================================================

'   Procedures:
'       <none>

'   Parameter:
'        Mean:  E[X] = expected value of  X

'   Local variables:
'       Index:  loop counter, second subscript of array  P( , ) .
'           n:  sample size (natural number power of 2)
'           i:  loop counter
'       Total:  loop counter
'        Prob:  temporary variable for  P(i, 1)
'         Sum:  Sum_i P(i, 1) (must = 1) ; <later:  sum of backward diagonal>
'         Min:  initial value of loop for calculating  P(Total,Index) .
'         Max:  final value of loop for calculating  P(Total,Index) .

DEFINT I, N
DIM Min AS INTEGER, Max AS INTEGER
CLS
PRINT TAB(20); "Illustration for the Central Limit Theorem"
PRINT "A six-sided die is loaded in accordance with probabilities that you"
PRINT "specify.   The probability mass function for the mean  Xbar  of a random"
PRINT "sample of size  n  will then be produced for values of  n  which will"
PRINT "double each time you press any key.   The mean value  E[X] is highlighted."
PRINT
PRINT "Working out probabilities for 1 die"

DO
    FOR i = 1 TO 6                                          ' Single roll.
        LOCATE 8 + i, 1
        PRINT SPC(75); " "
        LOCATE 8 + i, 1
        PRINT USING "Enter the probability of a score of # "; i;
        INPUT "on the die: ", Prob
       
        DO WHILE Prob < 0 OR Prob > 1            ' Reject invalid probability.
            LOCATE 9 + i, 1
            PRINT "Invalid value.   Must have  0 <= prob <= 1 .   ";
            INPUT "Try again: ", Prob
        LOOP
       
        LOCATE 9 + i, 1
        PRINT SPC(75); " "                       ' Erase any error message.
        LOCATE 8 + i, 51
        PRINT Prob
        LET P(i, 1) = Prob                       ' Record valid probability.
    NEXT i

    LET Sum = 0
   
    FOR i = 1 TO 6
        LET Sum = Sum + P(i, 1)
    NEXT i
   
    IF ABS(Sum - 1) >= .00001 THEN     ' Check that probabilities add up to 1.
        PRINT
        PRINT "Probabilities don't add up to 1!    Redo from start."
    ELSE
        LOCATE 15, 1                   ' Wipe out any error message.
        PRINT SPC(75); " "
        PRINT SPC(75); " "
        PRINT SPC(75); " "
        PRINT SPC(75); " "
    END IF

LOOP UNTIL ABS(Sum - 1) < .00001

LET Mean = 0                           ' Calculate the mean value,  E[X] .

FOR i = 1 TO 6
    LET Mean = Mean + i * P(i, 1)
NEXT i

'--------------------------------------   Now work out  P(Xbar = x)  for
'                                                   various sample sizes  n .

LET n = 1

FOR Index = 1 TO 6
    LET n = 2 * n                      ' i.e. for n = 2, 4, 8, 16, 32 and 64

    LOCATE 19, 1
    PRINT USING "Working out probabilities for ## dice"; n

    FOR Total = n TO 3 * n + n / 2           ' One half of table.
        LOCATE 21, 1
        PRINT USING "Evaluating  P(dice total = ###)"; Total
        LET Sum = 0
        LET Min = n / 2
        LET Max = Total - Min

        FOR i = Min TO Max                   ' Sum along backward diagonal.
            LET Sum = Sum + P(i, Index) * P(Total - i, Index)
        NEXT i

        LET P(Total, Index + 1) = Sum
    NEXT Total

    FOR Total = 3 * n + n / 2 + 1 TO 6 * n    ' Other half of table.
        LOCATE 21, 1
        PRINT USING "Evaluating  P(dice total = ###)"; Total
        LET Sum = 0
        LET Max = 3 * n
        LET Min = Total - Max
       
        FOR i = Min TO Max                   ' Sum along backward diagonal.
            LET Sum = Sum + P(i, Index) * P(Total - i, Index)
        NEXT i

        LET P(Total, Index + 1) = Sum
    NEXT Total

NEXT Index
       
END SUB

DEFSNG I, N
FUNCTION f (x, n AS INTEGER)
'============================================================================
'   Finds  P(Xbar = x | sample size = n).
'============================================================================

'   Local variables:
'       T:  total of all  n  rolls;   Xbar = T / n .
'    area:  area of histogram bar = P(Xbar = x);  (f = height, width = 1/n)

'   Procedures:
'       <none>

DIM T AS INTEGER

IF n = 1 THEN
    LET f = P(x, 1)
ELSE
    LET T = CINT(n * x)                     ' Guard against roundoff errors.
    LET area = P(T, 1 + CINT(LOG(n) / LOG(2)))  ' n = 2^(index-1)
    LET f = area * n
END IF

END FUNCTION

SUB Histogram (n AS INTEGER, Mean)
'============================================================================
'  Draw an histogram for the sample mean  X_n
'============================================================================

'   Local Variables:
'       Class:  loop counter for current histogram interval.
'           x:  corresponding value of sample mean  Xbar .
'    Response:  the key pressed by the user.

'   Subprograms / functions:
'       IntroScreen:    (re)print welcome message to the user.
'            f(x,n):    find  P(Xbar = x)  with the sample size  n .

DIM Class AS INTEGER, Response AS STRING

CALL IntroScreen
LOCATE 7, 70
PRINT "n ="; n

FOR Class = 0 TO 5 * n
    LET x = 1 + Class / n
    LINE ((64 * x + 67 - 32 / n), 180)-STEP(64 / n, -80 * f(x, n)), , BF
                                        ' Draw histogram bar        --^
NEXT Class

LINE ((64 * Mean + 66), 180)-STEP(2, -150), 0, BF   ' Black out the mean.
LOCATE 24, 9 + CINT(8 * Mean)
PRINT CHR$(230);                        ' Print mu on the horizontal axis.
LOCATE 2, 11
PRINT "Press <space> for next n; "; CHR$(34); "-"; CHR$(34); " to go back; ";
PRINT CHR$(34); "X"; CHR$(34); " to exit"


    DO                                  ' Pause until the user hits any key.
        LET Response = UCASE$(INKEY$)
    LOOP UNTIL Response <> ""

    SELECT CASE ASC(Response)
        CASE 13, 32, 70, 78             ' <CR>, < >, F,N,f,n  all = next n.
            IF n < 64 THEN LET n = 2 * n ' (Can't go to n > 64)
        CASE 45, 66, 80                 ' - , B,P,b,P  =  previous n
            IF n > 1 THEN LET n = n / 2 ' (Can't go to n < 1)
        CASE 81, 88                     ' Q,q,X,x = quit program
            LET n = 100                 ' Triggers exit from program
    END SELECT                          ' Else don't change n (repeats
                                        '  the current graph).


IF n < 100 THEN CLS                     ' Clear screen for next graph.

END SUB

SUB IntroScreen
'============================================================================
'  On new graphics screen print a reminder of the chosen probabilities for
'  one roll of the die;  and set up the axes for the  p.m.f.  histogram.
'============================================================================

'   Local variables:
'          i:   loop counter.
'       Tick:   loop counter to place tick marks on each axis.

DIM i AS INTEGER, Tick AS INTEGER

SCREEN 2
PALETTE 0, 15
PALETTE 1, 0
CLS
PRINT TAB(15); "Illustration for the Central Limit Theorem"

FOR i = 1 TO 6                             ' Echo probabilities selected by
    LOCATE 3 + i, 1                                  '              the user.
    PRINT USING "P(#) ="; i;
    PRINT P(i, 1)
NEXT i

LINE (103, 70)-(103, 180)                            ' vertical axis
LINE (103, 180)-(500, 180)                           ' horizontal axis

FOR Tick = 1 TO 6
    LINE (3 + (Tick + 1) * 64, 178)-STEP(0, 4)       ' tick marks <--+-->
    LOCATE 24, 8 + 8 * Tick
    PRINT Tick;
    LINE (100, 180 - 16 * Tick)-STEP(6, 0)         ' tick marks vertical axis
    LOCATE 23 - 2 * Tick, 9
    PRINT USING "#.#"; Tick / 5
NEXT Tick

LOCATE 22, 65
PRINT "_"                                            ' horizontal axis label.
LOCATE 23, 65
PRINT "x";

END SUB