// This program draws the normal curve or "bell curve" used in statistics. // You can pass in alternate number of points (default is 1000) to generate // curves for different values. use statistics.frink plotNormal[mean, sigma, steps, g is graphics] := { low = 1/steps // Use rational numbers so that the exactly high = 1-low // right number of points is plotted. println["low is \$low"] minSigma = inversePhi[low, 8] println["minsigma is \$minSigma"] maxSigma = inversePhi[high, 8] println["maxsigma is \$maxSigma"] vscale = 8 sigma^2 // Found experimentally to look good. ceilingH = normalDensity[mean + sigma * maxSigma, mean, sigma] scaledCeilingH = ceilingH * vscale r = scaledCeilingH println["ceiling H is \$ceilingH"] println["scaled ceiling is \$scaledCeilingH"] g.color[0.5,0.5,0.5] g.line[mean + (minSigma * sigma), 0, mean + (maxSigma * sigma), 0] // Draw SD lines highestSD = floor[maxSigma] for s = -highestSD to highestSD { x = mean + (sigma * s) y = -normalDensity[x, mean, sigma] * vscale g.line[x, 0, x, y] } width = maxSigma - minSigma // This polyline is the normal curve. c = new polyline for s=minSigma to maxSigma+0.001 step (width/100) { x = mean + (sigma * s) y = -normalDensity[x, mean, sigma] * vscale c.addPoint[x,y] } g.add[c] g.color[0,0,0] wheel = r/2 first = true points = 0 for phi = high to low step ((low-high)/(steps-1)) { // s = now[] x = inversePhi[phi,100,15] n = normalDensity[x, mean, sigma] do { wheel = (wheel + 0.618034) mod 1 } while wheel > n h = wheel if first { g.color[1,0,0] // Draw the "you" circle in red. g.fillEllipseCenter[x, -1/2 r, r, r] g.color[0,0,0] g.font["SansSerif", 4] g.text["You are here.", x, 7] g.line[x, 5, x, 1] // Arrow body // Arrowhead p=new filledPolygon p.addPoint[x,.65] p.addPoint[x+0.3,2.5] p.addPoint[x-0.3,2.5] g.add[p] first = false } else g.fillEllipseCenter[x, -h*vscale, r, r] points = points+1 // e = now[] // println["point \$points, time is " + format[e-s,"ms",3]] } println["\$points points plotted."] } g = new graphics points = 1000 // You can pass in a number of points as the sole argument. if length[ARGS] > 0 points = eval[ARGS@0] plotNormal[100, 15, points, g] g.show[] g.write["normal\$points.svg", 1024, undef] g.write["normal\$points.png", 2000, undef] g.write["normal\$points.html", 800, undef] //g.print[]