normalStats(100)
print
normalStats(1000)
print
normalStats(10000)
print
normalStats(100000)
end

sub randomNormal()
    local u1, u2
	u1 = ran()
    u2 = ran()
    return cos(2 * pi * u1) * sqrt(-2 * log(u2))
end sub

sub normalStats(n)
    local r(n), h(11)
    local i, sum, hSum, mean, sd, scale, stars
	sum = 0
    hSum = 0

    // Generate samples: mean 0.5, std 0.25
    for i = 1 to n
        r(i) = 0.5 + randomNormal() / 4
        sum = sum + r(i)

        if r(i) < 0 then
            h(0) = h(0) + 1
        elsif r(i) >= 1 then
            h(11) = h(11) + 1
        else
            h(int(r(i)*10)+1) = h(int(r(i)*10)+1) + 1
        fi
    next i

    for i = 0 to 11
        hSum = hSum + h(i)
    next i

    mean = sum / n

    // Standard deviation
    sum = 0
    for i = 1 to n
        sum = sum + (r(i) - mean)^2
    next i
    sd = sqrt(sum / n)

    print "Sample size ", n
    print
    print "  Mean ", mean using("#.######");
    print "  SD ", sd using("#.######")
    print

    // Histogram (scale for large samples)
    scale = 1
    if n > 300  scale = 300 / n

    for i = 0 to 11
        if i = 0 then
            print "< 0.00  : ", h(0) using("#####"), " ";
        elsif i = 11 then
            print ">=1.00  : ", h(1) using("#####");
        else
            //print using("  #.##  : ##### ", (i-1)/10); h(i);
			print "  ", (i-1)/10 using("#.##"), "  : ", h(i) using("#####"), " ";
        fi

        stars = int(h(i) * scale + 0.5)
        print string$(stars, "*")
    next i
end sub