1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
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);
|
