' Random distributions: uniform, normal, exponential, Bernoulli, Poisson
#COMPILE EXE
#DIM ALL

FUNCTION UniformRand(BYVAL min AS DOUBLE, BYVAL max AS DOUBLE) AS DOUBLE
    FUNCTION = min + (max - min) * RND
END FUNCTION

FUNCTION NormalRand(BYVAL mean AS DOUBLE, BYVAL stddev AS DOUBLE) AS DOUBLE
    LOCAL u1  AS DOUBLE
    LOCAL u2  AS DOUBLE
    LOCAL r  AS DOUBLE
    LOCAL theta  AS DOUBLE
    u1 = RND
    IF u1 < 0.000000000000001 THEN u1 = 0.000000000000001
    u2 = RND
    r = SQR(0.0 - 2.0 * 1)
    theta = 6.283185307179586 * u2
    FUNCTION = mean + stddev * r * COS(theta)
END FUNCTION

FUNCTION ExponentialRand(BYVAL lambda AS DOUBLE) AS DOUBLE
    LOCAL u AS DOUBLE
    u = RND
    IF u < 0.000000000000001 THEN u = 0.000000000000001
    FUNCTION = (0.0 - 1) / lambda
END FUNCTION

FUNCTION BernoulliRand(BYVAL p AS DOUBLE) AS LONG
    IF RND < p THEN FUNCTION = 1 ELSE FUNCTION = 0
END FUNCTION

FUNCTION PoissonRand(BYVAL lambda AS DOUBLE) AS LONG
    LOCAL L  AS DOUBLE
    LOCAL p  AS DOUBLE
    LOCAL u  AS DOUBLE
    LOCAL k AS LONG
    L = (2.718281828459045 ^ (-lambda))
    k = 0
    p = 1.0
    DO
        k = k + 1
        u = RND
        IF u < 0.000000000000001 THEN u = 0.000000000000001
        p = p * u
        IF p < L THEN EXIT LOOP
    LOOP
    FUNCTION = k - 1
END FUNCTION

FUNCTION SampleMean(BYVAL n AS LONG) AS DOUBLE
    LOCAL i AS LONG
    LOCAL sum AS DOUBLE
    sum = 0.0
    FOR i = 1 TO n
        sum = sum + UniformRand(0.0, 1.0)
    NEXT i
    FUNCTION = sum / n
END FUNCTION

FUNCTION PBMAIN() AS LONG
    LOCAL i  AS LONG
    LOCAL n  AS LONG
    LOCAL ones  AS LONG
    PRINT "=== Random Distributions ==="
    PRINT
    RANDOMIZE TIMER
    n = 10
    PRINT "Uniform[0,1]:"
    FOR i = 1 TO 5
        PRINT "  " + STR$(UniformRand(0.0, 1.0))
    NEXT i
    PRINT
    PRINT "Normal(0):"
    FOR i = 1 TO 5
        PRINT "  " + STR$(NormalRand(0.0, 1.0))
    NEXT i
    PRINT
    PRINT "Exponential(1):"
    FOR i = 1 TO 5
        PRINT "  " + STR$(ExponentialRand(1.0))
    NEXT i
    PRINT
    ones = 0
    FOR i = 1 TO n
        IF BernoulliRand(0.5) THEN ones = ones + 1
    NEXT i
    PRINT "Bernoulli(0.5): " + STR$(ones) + " ones / " + STR$(n)
    PRINT
    PRINT "Poisson(3):"
    FOR i = 1 TO 5
        PRINT "  " + STR$(PoissonRand(3.0))
    NEXT i
    PRINT
    PRINT "Uniform(10,20): " + STR$(UniformRand(10.0, 20.0))
    PRINT "Normal(100,15): " + STR$(NormalRand(100.0, 15.0))
    PRINT "Sample mean(100): " + STR$(SampleMean(100))
    PRINT
    PRINT "Distribution tests complete."
    FUNCTION = 0
END FUNCTION
