DECLARE FUNCTION randomNormal! ()
DECLARE SUB normalStats (sampleSize)

CONST pi = 3.141592653589793#
RANDOMIZE TIMER

CLS
CALL normalStats(100)
PRINT
CALL normalStats(1000)
PRINT
CALL normalStats(10000)
PRINT
CALL normalStats(100000)
PRINT "Press any key to quit"
END

SUB normalStats (sampleSize)
    IF sampleSize < 1 THEN EXIT SUB

    DIM r(1 TO sampleSize) AS SINGLE
    DIM h(-1 TO 10) AS INTEGER
    DIM sum AS DOUBLE: sum = 0!
    DIM hSum AS INTEGER: hSum = 0
    DIM i AS INTEGER

    ' Generate samples: mean 0.5, std 0.25
    FOR i = 1 TO sampleSize
        r(i) = .5 + randomNormal! / 4!
        sum = sum + r(i)

        IF r(i) < 0! THEN
            h(-1) = h(-1) + 1
        ELSEIF r(i) >= 1! THEN
            h(10) = h(10) + 1
        ELSE
            h(INT(r(i) * 10)) = h(INT(r(i) * 10)) + 1
        END IF
    NEXT i

    FOR i = -1 TO 10
        hSum = hSum + h(i)
    NEXT i

    ' Adjust if necessary (rounding error)
    DIM adj AS INTEGER: adj = sampleSize - hSum
    IF adj <> 0 THEN
        FOR i = -1 TO 10
            h(i) = h(i) + adj
            IF h(i) >= 0 THEN EXIT FOR
            h(i) = h(i) - adj
        NEXT i
    END IF

    DIM mean AS DOUBLE: mean = sum / sampleSize

    ' Standard deviation
    sum = 0!
    FOR i = 1 TO sampleSize
        sum = sum + (r(i) - mean) ^ 2!
    NEXT i
    DIM sd AS DOUBLE: sd = SQR(sum / sampleSize)

    PRINT "Sample size"; sampleSize
    PRINT
    PRINT USING "  Mean #.######"; mean;
    PRINT USING "  SD #.######"; sd
    PRINT

    ' Histogram (scale large samples)
    DIM numStars AS INTEGER
    DIM scale AS DOUBLE: scale = 1!
    IF sampleSize > 300 THEN scale = 300 / sampleSize

    FOR i = -1 TO 10
        IF i = -1 THEN
            PRINT "< 0.00  : ";
        ELSEIF i = 10 THEN
            PRINT ">=1.00  : ";
        ELSE
            PRINT USING "  #.##  : "; i / 10;
        END IF

        PRINT USING "##### "; h(i);
        numStars = INT(h(i) * scale + .5)
        PRINT STRING$(numStars, "*")
    NEXT i
END SUB

FUNCTION randomNormal!
    randomNormal! = COS(2 * pi * RND) * SQR(-2 * LOG(RND))
END FUNCTION