diff options
Diffstat (limited to 'Build/source/utils/asymptote/examples/linearregression.asy')
-rw-r--r-- | Build/source/utils/asymptote/examples/linearregression.asy | 64 |
1 files changed, 64 insertions, 0 deletions
diff --git a/Build/source/utils/asymptote/examples/linearregression.asy b/Build/source/utils/asymptote/examples/linearregression.asy new file mode 100644 index 00000000000..a1cc6ff6bcb --- /dev/null +++ b/Build/source/utils/asymptote/examples/linearregression.asy @@ -0,0 +1,64 @@ +import graph3; +import math; // for the leastsquares routine + +Billboard.targetsize = true; // Perspective should not affect the labels. +currentprojection = perspective(60 * (5, 2, 3)); + +file duncan = input("linearregression.dat"); + +string headers = duncan; + +real[][] independentvars; +real[] dependentvars; + +while (!eof(duncan)) { + string line = duncan; + string[] entries = split(line); + if (entries.length < 5) continue; + string type = entries[1]; + real income = (real)(entries[2]); + real education = (real)(entries[3]); + real prestige = (real)(entries[4]); + + // include 1.0 for the residue + independentvars.push(new real[] {income, education, 1.0}); + dependentvars.push(prestige); +} + +real[] coeffs = leastsquares(independentvars, dependentvars, warn=false); +if (coeffs.length == 0) { + abort("Unable to find regression: independent variables are " + + "linearly dependent."); +} + +real f(pair xy) { + return coeffs[0] * xy.x // income + + coeffs[1] * xy.y // education + + coeffs[2]; // residue +} + +real xmin = infinity, xmax = -infinity, ymin = infinity, ymax = -infinity; +for (real[] row : independentvars) { + if (row[0] < xmin) xmin = row[0]; + if (row[0] > xmax) xmax = row[0]; + if (row[1] < ymin) ymin = row[1]; + if (row[1] > ymax) ymax = row[1]; +} + +// Draw the plane +draw(surface(f, (xmin, ymin), (xmax, ymax)), + surfacepen=emissive(blue + opacity(0.6)), + meshpen = blue); + +for (int ii = 0; ii < independentvars.length; ++ii) { + triple pt = (independentvars[ii][0], independentvars[ii][1], + dependentvars[ii]); + draw(shift(pt) * unitsphere, material(yellow, emissivepen=0.2*yellow)); + real z = f((pt.x, pt.y)); + if (pt.z > z) draw (pt -- (pt.x, pt.y, z), green); + else draw(pt -- (pt.x, pt.y, z), red); +} + +xaxis3("income", Bounds(Min, Min), InTicks); +yaxis3("education", Bounds(Min, Min), InTicks); +zaxis3("prestige", Bounds(Min, Min), InTicks); |