! Fisher's Discriminant analysis (FisherDA.lng);
! There are two categories of observations, type 0 and type 1,
which we wish to distinguish based on explanatory variables;
! Key idea: Find a scoring formula:
score( i) = beta0 + beta1*xdat( k, 1) + beta2*xdat( i, 2) + ...
so to extent possible:
score( i) < 0 if observation i is of type 0,
score( i) > 0 if observation i is of type 1;
! Reference: Fisher, R.A. (1936), “The Use of Multiple Measurements in Taxonomical Problems,”
Annals of Eugenics, vol. 7, pp. 179-188.;
!Keywords: Binary choice, Categorical variables, Discriminant analysis, Fisher, Random forest;
SETS:
var: BETA,
mu0, MU1, mutot ;
vardep( VAR); ! The dependent variable;
varexp( VAR); ! The set of explanatory variables;
varpltx( var); ! Variable on X axis in plot/chart/graph;
varplty( var); ! Variable on Y axis in plot/chart/graph;
varwgt( var); ! The weight vector column (if used);
OBS: ERR, SCORE;
oxv( OBS, var): xdat;
vxv( var, var): covar, covarinv;
! Two subsets, based on category of the observation;
OBS0( OBS) : xp0, yp0;
OBS1( OBS) : xp1, yp1;
ENDSETS
DATA:
! Various data sets;
! Case: a small circle, Type 0, inside large circle, Type 1;
!CaCircle var = X1 X2 X1SQ X2SQ type wgtof;
!CaCircle vardep = type;
!CaCircle !CaCircle varexp = X1 X2 X1SQ X2SQ; !Can if include squared terms!CaCircle varpltx = X1; !CaCircle varplty = X2; !CaCircle varwgt = wgtof; !CaCircle xdat =
22.98 21.75 528.30 473.06 0 1
28.27 24.90 799.05 620.01 1 1
21.85 22.88 477.42 523.71 0 1
24.90 28.34 620.01 803.19 1 1
20.30 23.30 412.09 542.89 0 1
20.30 29.60 412.09 876.16 1 1
18.75 22.88 351.56 523.71 0 1
15.70 28.34 246.49 803.19 1 1
17.62 21.75 310.30 473.06 0 1
12.33 24.90 152.09 620.01 1 1
17.20 20.20 295.84 408.04 0 1
11.10 20.20 123.21 408.04 1 1
17.62 18.65 310.30 347.82 0 1
12.33 15.50 152.09 240.25 1 1
18.75 17.52 351.56 306.79 0 1
15.70 12.06 246.49 145.43 1 1
20.30 17.10 412.09 292.41 0 1
20.30 10.80 412.09 116.64 1 1
21.85 17.52 477.42 306.79 0 1
24.90 12.06 620.01 145.43 1 1
22.98 18.65 528.30 347.82 0 1
28.27 15.50 799.05 240.25 1 1
23.40 20.20 547.56 408.04 0 1
29.50 20.20 870.25 408.04 1 1;
! Case 1 with Unequal class sizes;
!CaUneq01 var = value type wgtof; !CaUneq01 vardep = type; !CaUneq01 varexp = value; !CaUneq01 varpltx = value; !CaUneq01 varplty = value;
!CaUneq01 xdat =
1 1 1
2 0 1
3 0 1
7 0 1
8 0 1
9 0 1;
! Case 2 Unequal class sizes;
!CaUneq02 var = value type wgtof;
!CaUneq02 vardep = type;
!CaUneq02 varexp = value; !CaUneq02 varpltx = value;
!CaUneq02 varplty = value;
!CaUneq02 varwgt = wgtof;
!CaUneq02 xdat =
1 0 1
2 0 1
3 1 1
7 0 1
8 0 1
9 0 1;
! An inseparable case;
!CaInsep var = value type wgtof;
!CaInsep vardep = type;
!CaInsep varexp = value;
!CaInsep varpltx = value; !CaInsep varplty = wgtof; !CaInsep varwgt = wgtof;
!CaInsep alpha = 0; !CaInsep xdat =
1 0 1
2 0 1
2 1 1
3 1 1;
! A separable case;
!CaFDA var = X1 X2 type wgtof;
!CaFDA vardep = type;
!CaFDA varexp = X1 X2; !CaFDA varpltx = X1; !CaFDA varplty = X2; !CaFDA varwgt = wgtof;
!CaFDA alpha = 0; !CaFDA xdat =
4 1 0 1
2 4 0 1
2 3 0 1
3 6 0 1
4 4 0 1
9 10 1 1
6 8 1 1
9 5 1 1
8 7 1 1
10 8 1 1;
! A separable case, csd.uwo.ca, should get (beta0, beta1) ~ (-0.79, 0.89);
!CaFDA01 var = X1 X2 type wgtof;
!CaFDA01 vardep = type; !CaFDA01 varexp = X1 X2; !CaFDA01 varpltx = X1; !CaFDA01 varplty = X2; !CaFDA01 varwgt = wgtof;
!CaFDA01 alpha = 0; !CaFDA01 xdat =
1 2 0 1
2 3 0 1
3 3 0 1
4 5 0 1
5 5 0 1
1 0 1 1
2 1 1 1
3 1 1 1
3 2 1 1
5 3 1 1
6 5 1 1;
! Unequal class sizes, case 3;
!CaUneq03 var = value type wgtof;
!CaUneq03 vardep = type;
!CaUneq03 varexp = value; !CaUneq03 varpltx = value; !CaUneq03 varplty = wgtof; !CaUneq03 varwgt = wgtof;
!CaUneq03 xdat =
1 0 0.6
2 0 0.6
3 1 3
7 0 0.6
8 0 0.6
9 0 0.6;
! Illustrate separability and parsimony;
!CaSePar var = X1 X2 type wgtof;
!CaSePar vardep = type;
!CaSePar varexp = X1 X2; !CaSePar varpltx = X1; !CaSePar varplty = X2; !CaSePar varwgt = wgtof;
!CaSePar xdat =
1 8 0 1
2 6 0 1
3 9 0 1
4 7 0 1
6 3 1 1
7 1 1 1
8 4 1 1
9 2 1 1;
! The 66 variable set;
! OBS = 1..66;
!Ca66 var =
OBNO TEST1 TEST2 TEST3 TEST4 TEST5 TYPE;
!Ca66 vardep = type;
!Ca66 varexp=
TEST1 TEST2 TEST3 TEST4 TEST5;
!Ca66 varpltx = TEST1;
!Ca66 varplty = TEST2;
! Variable on X axis in plot/chart/graph;
!Ca66 xdat =
1 346.7000122 247.1999969 220.5 364.1000061 311.7000122 0
2 334 313.2999878 306.5 330.8999939 311.1000061 0
3 248.3999939 189.1999969 206.8000031 334.7000122 312.5 1
4 309 291.8999939 281.2000122 346.2000122 311.1000061 1
5 328.8999939 306.2000122 259.3999939 336.3999939 310.8999939 0
6 252.8000031 248.8000031 253.8000031 321 311.7000122 1
7 313 289.7000122 292.6000061 318 311 0
8 304.8999939 115.5 284.2000122 316.5 310.5 1
9 327.8999939 330.7999878 305.7000122 332.6000061 311 1
10 315.3999939 203.8999939 287.1000061 333.7999878 311.5 1
11 333 270.6000061 274.2999878 379.1000061 311.2000122 0
12 242.3999939 145.8999939 292.2999878 318.7000122 311.2999878 0
13 124.9000015 1.1 244.1999969 345.7000122 310.7999878 0
14 323.5 317.2000122 287.3999939 406.1000061 312 1
15 304.2999878 191.6999969 275.7999878 331.7000122 311.5 0
16 382.3999939 124.0999985 30 322.5 316.7000122 1
17 327 275.3999939 290.6000061 345.5 313.3999939 1
18 278.7999878 282.1000061 316.2999878 317 311.2999878 0
19 324.1000061 261.7999878 316.7999878 326.6000061 311.6000061 0
20 249.3999939 260.7999878 292.7999878 317.2000122 310.2999878 1
21 336.2000122 290.7999878 273.2999878 400.3999939 310.7999878 1
22 317 291.8999939 303.5 326.5 310.8999939 0
23 256.8999939 212 282 336.6000061 311.7000122 1
24 292.7999878 181 295.7999878 577.9000244 311.2999878 0
25 342.7000122 306 294.2000122 487.3999939 312.1000061 0
26 336.7000122 301.2999878 273.7000122 342.5 312.7999878 0
27 302.2999878 250.8000031 297.2000122 331.2999878 312.1000061 0
28 328 296.8999939 292.3999939 324.6000061 310.8999939 1
29 312 272 311.6000061 317.7000122 311.2000122 1
30 274.7000122 252.1000061 310.7000122 323.7000122 310.7999878 0
31 315.1000061 301.2000122 300.8999939 410.8999939 310.8999939 0
32 310 245.3000031 306 310.7000122 310.1000061 0
33 335.2000122 298.6000061 314.7999878 317 310.8999939 1
34 345.2000122 353 326.3999939 409.1000061 311.2999878 1
35 348.7999878 357 326 436.5 311.8999939 0
36 324 306.7000122 314 401.7000122 312.7000122 0
37 365.1000061 345 330.7999878 382.2999878 311.8999939 1
38 369.2999878 356.7000122 322.6000061 1034.1 310.8999939 1
39 343.6000061 330.7999878 322.5 462.7999878 312.3999939 1
40 362.7999878 343 333.6000061 785.9000244 311.5 1
41 355.6000061 336.1000061 320.3999939 597.9000244 312.1000061 0
42 357.3999939 378.6000061 323.7999878 891.2999878 311.6000061 1
43 350 347.2999878 343.3999939 538.7999878 313.5 1
44 379 369 333.1000061 716 315.5 1
45 344.2000122 359.6000061 333.7999878 436.6000061 311.8999939 1
46 357 322.5 317 363.3999939 311.7999878 0
47 325.3999939 347.2999878 344.1000061 880.0999756 311.5 1
48 366.8999939 345.2999878 314.2000122 550.2999878 310.8999939 0
49 353.7999878 359.5 335.1000061 425 312.6000061 0
50 330.7000122 328.1000061 323.5 373.1000061 314 1
51 343.7999878 341.3999939 325.7000122 454.7999878 311.8999939 0
52 345.2999878 331.5 295.6000061 400 311 1
53 334.3999939 318.5 315.7999878 459.1000061 311.5 0
54 358.8999939 350.6000061 315.7999878 392 311.7999878 0
55 359.8999939 344.6000061 336.3999939 620 311.7999878 0
56 364.7999878 329.8999939 336.7000122 549.9000244 312.2999878 1
57 349 327.3999939 322.6000061 370.5 311.2999878 1
58 363 364.7000122 324.6000061 1081.7 311.7000122 1
59 330.1000061 363.5 330.6000061 617.5 311.1000061 1
60 363.7000122 345.8999939 336.3999939 599.5 312 0
61 356.1000061 349.3999939 340.5 1010 311.8999939 1
62 358.2999878 363.1000061 317.1000061 474.3999939 311.8999939 0
63 356.7000122 349.7999878 323.7999878 539.0999756 311.2000122 0
64 370.2999878 369.5 317 536.5999756 312 1
65 327.8999939 326.2999878 330.3999939 415.6000061 311 0
66 334.7000122 331.7000122 302.2000122 428.6000061 311.6000061 0;
!Fisher's Iris Data set. setosa is easily distinguished from the other 2.
Observations 1-50 are setosa, 51-100 are versicolor, 101-150 are virginica.
!CaIris var = Order Sepallen SepalWid PetalLen PetalWid Species wgtof;
!CaIris vardep = Species; !CaIris varexp = Sepallen SepalWid PetalLen PetalWid;
!CaIris varpltx = Sepallen ; !CaIris varplty = SepalWid ; !CaIris varwgt = wgtof; !Data Sepal Sepal Petal Petal
Order length width length width Species Wgt;
!CaIris xdat =
1 5.1 3.5 1.4 0.2 0 1
2 4.9 3.0 1.4 0.2 0 1
3 4.7 3.2 1.3 0.2 0 1
4 4.6 3.1 1.5 0.2 0 1
5 5.0 3.6 1.4 0.3 0 1
6 5.4 3.9 1.7 0.4 0 1
7 4.6 3.4 1.4 0.3 0 1
8 5.0 3.4 1.5 0.2 0 1
9 4.4 2.9 1.4 0.2 0 1
10 4.9 3.1 1.5 0.1 0 1
11 5.4 3.7 1.5 0.2 0 1
12 4.8 3.4 1.6 0.2 0 1
13 4.8 3.0 1.4 0.1 0 1
14 4.3 3.0 1.1 0.1 0 1
15 5.8 4.0 1.2 0.2 0 1
16 5.7 4.4 1.5 0.4 0 1
17 5.4 3.9 1.3 0.4 0 1
18 5.1 3.5 1.4 0.3 0 1
19 5.7 3.8 1.7 0.3 0 1
20 5.1 3.8 1.5 0.3 0 1
21 5.4 3.4 1.7 0.2 0 1
22 5.1 3.7 1.5 0.4 0 1
23 4.6 3.6 1.0 0.2 0 1
24 5.1 3.3 1.7 0.5 0 1
25 4.8 3.4 1.9 0.2 0 1
26 5.0 3.0 1.6 0.2 0 1
27 5.0 3.4 1.6 0.4 0 1
28 5.2 3.5 1.5 0.2 0 1
29 5.2 3.4 1.4 0.2 0 1
30 4.7 3.2 1.6 0.2 0 1
31 4.8 3.1 1.6 0.2 0 1
32 5.4 3.4 1.5 0.4 0 1
33 5.2 4.1 1.5 0.1 0 1
34 5.5 4.2 1.4 0.2 0 1
35 4.9 3.1 1.5 0.2 0 1
36 5.0 3.2 1.2 0.2 0 1
37 5.5 3.5 1.3 0.2 0 1
38 4.9 3.6 1.4 0.1 0 1
39 4.4 3.0 1.3 0.2 0 1
40 5.1 3.4 1.5 0.2 0 1
41 5.0 3.5 1.3 0.3 0 1
42 4.5 2.3 1.3 0.3 0 1
43 4.4 3.2 1.3 0.2 0 1
44 5.0 3.5 1.6 0.6 0 1
45 5.1 3.8 1.9 0.4 0 1
46 4.8 3.0 1.4 0.3 0 1
47 5.1 3.8 1.6 0.2 0 1
48 4.6 3.2 1.4 0.2 0 1
49 5.3 3.7 1.5 0.2 0 1
50 5.0 3.3 1.4 0.2 0 1
51 7.0 3.2 4.7 1.4 1 1
52 6.4 3.2 4.5 1.5 1 1
53 6.9 3.1 4.9 1.5 1 1
54 5.5 2.3 4.0 1.3 1 1
55 6.5 2.8 4.6 1.5 1 1
56 5.7 2.8 4.5 1.3 1 1
57 6.3 3.3 4.7 1.6 1 1
58 4.9 2.4 3.3 1.0 1 1
59 6.6 2.9 4.6 1.3 1 1
60 5.2 2.7 3.9 1.4 1 1
61 5.0 2.0 3.5 1.0 1 1
62 5.9 3.0 4.2 1.5 1 1
63 6.0 2.2 4.0 1.0 1 1
64 6.1 2.9 4.7 1.4 1 1
65 5.6 2.9 3.6 1.3 1 1
66 6.7 3.1 4.4 1.4 1 1
67 5.6 3.0 4.5 1.5 1 1
68 5.8 2.7 4.1 1.0 1 1
69 6.2 2.2 4.5 1.5 1 1
70 5.6 2.5 3.9 1.1 1 1
71 5.9 3.2 4.8 1.8 1 1
72 6.1 2.8 4.0 1.3 1 1
73 6.3 2.5 4.9 1.5 1 1
74 6.1 2.8 4.7 1.2 1 1
75 6.4 2.9 4.3 1.3 1 1
76 6.6 3.0 4.4 1.4 1 1
77 6.8 2.8 4.8 1.4 1 1
78 6.7 3.0 5.0 1.7 1 1
79 6.0 2.9 4.5 1.5 1 1
80 5.7 2.6 3.5 1.0 1 1
81 5.5 2.4 3.8 1.1 1 1
82 5.5 2.4 3.7 1.0 1 1
83 5.8 2.7 3.9 1.2 1 1
84 6.0 2.7 5.1 1.6 1 1
85 5.4 3.0 4.5 1.5 1 1
86 6.0 3.4 4.5 1.6 1 1
87 6.7 3.1 4.7 1.5 1 1
88 6.3 2.3 4.4 1.3 1 1
89 5.6 3.0 4.1 1.3 1 1
90 5.5 2.5 4.0 1.3 1 1
91 5.5 2.6 4.4 1.2 1 1
92 6.1 3.0 4.6 1.4 1 1
93 5.8 2.6 4.0 1.2 1 1
94 5.0 2.3 3.3 1.0 1 1
95 5.6 2.7 4.2 1.3 1 1
96 5.7 3.0 4.2 1.2 1 1
97 5.7 2.9 4.2 1.3 1 1
98 6.2 2.9 4.3 1.3 1 1
99 5.1 2.5 3.0 1.1 1 1
100 5.7 2.8 4.1 1.3 1 1
101 6.3 3.3 6.0 2.5 1 1
102 5.8 2.7 5.1 1.9 1 1
103 7.1 3.0 5.9 2.1 1 1
104 6.3 2.9 5.6 1.8 1 1
105 6.5 3.0 5.8 2.2 1 1
106 7.6 3.0 6.6 2.1 1 1
107 4.9 2.5 4.5 1.7 1 1
108 7.3 2.9 6.3 1.8 1 1
109 6.7 2.5 5.8 1.8 1 1
110 7.2 3.6 6.1 2.5 1 1
111 6.5 3.2 5.1 2.0 1 1
112 6.4 2.7 5.3 1.9 1 1
113 6.8 3.0 5.5 2.1 1 1
114 5.7 2.5 5.0 2.0 1 1
115 5.8 2.8 5.1 2.4 1 1
116 6.4 3.2 5.3 2.3 1 1
117 6.5 3.0 5.5 1.8 1 1
118 7.7 3.8 6.7 2.2 1 1
119 7.7 2.6 6.9 2.3 1 1
120 6.0 2.2 5.0 1.5 1 1
121 6.9 3.2 5.7 2.3 1 1
122 5.6 2.8 4.9 2.0 1 1
123 7.7 2.8 6.7 2.0 1 1
124 6.3 2.7 4.9 1.8 1 1
125 6.7 3.3 5.7 2.1 1 1
126 7.2 3.2 6.0 1.8 1 1
127 6.2 2.8 4.8 1.8 1 1
128 6.1 3.0 4.9 1.8 1 1
129 6.4 2.8 5.6 2.1 1 1
130 7.2 3.0 5.8 1.6 1 1
131 7.4 2.8 6.1 1.9 1 1
132 7.9 3.8 6.4 2.0 1 1
133 6.4 2.8 5.6 2.2 1 1
134 6.3 2.8 5.1 1.5 1 1
135 6.1 2.6 5.6 1.4 1 1
136 7.7 3.0 6.1 2.3 1 1
137 6.3 3.4 5.6 2.4 1 1
138 6.4 3.1 5.5 1.8 1 1
139 6.0 3.0 4.8 1.8 1 1
140 6.9 3.1 5.4 2.1 1 1
141 6.7 3.1 5.6 2.4 1 1
142 6.9 3.1 5.1 2.3 1 1
143 5.8 2.7 5.1 1.9 1 1
144 6.8 3.2 5.9 2.3 1 1
145 6.7 3.3 5.7 2.5 1 1
146 6.7 3.0 5.2 2.3 1 1
147 6.3 2.5 5.0 1.9 1 1
148 6.5 3.0 5.2 2.0 1 1
149 6.2 3.4 5.4 2.3 1 1
150 5.9 3.0 5.1 1.8 1 1
;
! Genuine and counterfeit banknotes (100 Swiss Francs),
various measurements.
Banknotes BN1 to BN100 are genuine (Good=1),
all others are counterfeit (Good=0).
Dataset courtesy of H. Riedwyl, Bern, Switzerland;
!CaSwiss; var =
Length Left Right Bottom Top Diagonal Good wgtof ;
!CaSwiss; vardep = Good!CaSwiss; varexp = Length Left Right Bottom Top Diagonal!CaSwiss; varpltx = Right!CaSwiss; varplty = Diagonal!CaSwiss; alpha = 0.0001!CaSwiss; varwgt = wgtof
!CaSwiss; OBS, xdat =
BN1 214.8 131.0 131.1 9.0 9.7 141.0 1 1
BN2 214.6 129.7 129.7 8.1 9.5 141.7 1 1
BN3 214.8 129.7 129.7 8.7 9.6 142.2 1 1
BN4 214.8 129.7 129.6 7.5 10.4 142.0 1 1
BN5 215.0 129.6 129.7 10.4 7.7 141.8 1 1
BN6 215.7 130.8 130.5 9.0 10.1 141.4 1 1
BN7 215.5 129.5 129.7 7.9 9.6 141.6 1 1
BN8 214.5 129.6 129.2 7.2 10.7 141.7 1 1
BN9 214.9 129.4 129.7 8.2 11.0 141.9 1 1
BN10 215.2 130.4 130.3 9.2 10.0 140.7 1 1
BN11 215.3 130.4 130.3 7.9 11.7 141.8 1 1
BN12 215.1 129.5 129.6 7.7 10.5 142.2 1 1
BN13 215.2 130.8 129.6 7.9 10.8 141.4 1 1
BN14 214.7 129.7 129.7 7.7 10.9 141.7 1 1
BN15 215.1 129.9 129.7 7.7 10.8 141.8 1 1
BN16 214.5 129.8 129.8 9.3 8.5 141.6 1 1
BN17 214.6 129.9 130.1 8.2 9.8 141.7 1 1
BN18 215.0 129.9 129.7 9.0 9.0 141.9 1 1
BN19 215.2 129.6 129.6 7.4 11.5 141.5 1 1
BN20 214.7 130.2 129.9 8.6 10.0 141.9 1 1
BN21 215.0 129.9 129.3 8.4 10.0 141.4 1 1
BN22 215.6 130.5 130.0 8.1 10.3 141.6 1 1
BN23 215.3 130.6 130.0 8.4 10.8 141.5 1 1
BN24 215.7 130.2 130.0 8.7 10.0 141.6 1 1
BN25 215.1 129.7 129.9 7.4 10.8 141.1 1 1
BN26 215.3 130.4 130.4 8.0 11.0 142.3 1 1
BN27 215.5 130.2 130.1 8.9 9.8 142.4 1 1
BN28 215.1 130.3 130.3 9.8 9.5 141.9 1 1
BN29 215.1 130.0 130.0 7.4 10.5 141.8 1 1
BN30 214.8 129.7 129.3 8.3 9.0 142.0 1 1
BN31 215.2 130.1 129.8 7.9 10.7 141.8 1 1
BN32 214.8 129.7 129.7 8.6 9.1 142.3 1 1
BN33 215.0 130.0 129.6 7.7 10.5 140.7 1 1
BN34 215.6 130.4 130.1 8.4 10.3 141.0 1 1
BN35 215.9 130.4 130.0 8.9 10.6 141.4 1 1
BN36 214.6 130.2 130.2 9.4 9.7 141.8 1 1
BN37 215.5 130.3 130.0 8.4 9.7 141.8 1 1
BN38 215.3 129.9 129.4 7.9 10.0 142.0 1 1
BN39 215.3 130.3 130.1 8.5 9.3 142.1 1 1
BN40 213.9 130.3 129.0 8.1 9.7 141.3 1 1
BN41 214.4 129.8 129.2 8.9 9.4 142.3 1 1
BN42 214.8 130.1 129.6 8.8 9.9 140.9 1 1
BN43 214.9 129.6 129.4 9.3 9.0 141.7 1 1
BN44 214.9 130.4 129.7 9.0 9.8 140.9 1 1
BN45 214.8 129.4 129.1 8.2 10.2 141.0 1 1
BN46 214.3 129.5 129.4 8.3 10.2 141.8 1 1
BN47 214.8 129.9 129.7 8.3 10.2 141.5 1 1
BN48 214.8 129.9 129.7 7.3 10.9 142.0 1 1
BN49 214.6 129.7 129.8 7.9 10.3 141.1 1 1
BN50 214.5 129.0 129.6 7.8 9.8 142.0 1 1
BN51 214.6 129.8 129.4 7.2 10.0 141.3 1 1
BN52 215.3 130.6 130.0 9.5 9.7 141.1 1 1
BN53 214.5 130.1 130.0 7.8 10.9 140.9 1 1
BN54 215.4 130.2 130.2 7.6 10.9 141.6 1 1
BN55 214.5 129.4 129.5 7.9 10.0 141.4 1 1
BN56 215.2 129.7 129.4 9.2 9.4 142.0 1 1
BN57 215.7 130.0 129.4 9.2 10.4 141.2 1 1
BN58 215.0 129.6 129.4 8.8 9.0 141.1 1 1
BN59 215.1 130.1 129.9 7.9 11.0 141.3 1 1
BN60 215.1 130.0 129.8 8.2 10.3 141.4 1 1
BN61 215.1 129.6 129.3 8.3 9.9 141.6 1 1
BN62 215.3 129.7 129.4 7.5 10.5 141.5 1 1
BN63 215.4 129.8 129.4 8.0 10.6 141.5 1 1
BN64 214.5 130.0 129.5 8.0 10.8 141.4 1 1
BN65 215.0 130.0 129.8 8.6 10.6 141.5 1 1
BN66 215.2 130.6 130.0 8.8 10.6 140.8 1 1
BN67 214.6 129.5 129.2 7.7 10.3 141.3 1 1
BN68 214.8 129.7 129.3 9.1 9.5 141.5 1 1
BN69 215.1 129.6 129.8 8.6 9.8 141.8 1 1
BN70 214.9 130.2 130.2 8.0 11.2 139.6 1 1
BN71 213.8 129.8 129.5 8.4 11.1 140.9 1 1
BN72 215.2 129.9 129.5 8.2 10.3 141.4 1 1
BN73 215.0 129.6 130.2 8.7 10.0 141.2 1 1
BN74 214.4 129.9 129.6 7.5 10.5 141.8 1 1
BN75 215.2 129.9 129.7 7.2 10.6 142.1 1 1
BN76 214.1 129.6 129.3 7.6 10.7 141.7 1 1
BN77 214.9 129.9 130.1 8.8 10.0 141.2 1 1
BN78 214.6 129.8 129.4 7.4 10.6 141.0 1 1
BN79 215.2 130.5 129.8 7.9 10.9 140.9 1 1
BN80 214.6 129.9 129.4 7.9 10.0 141.8 1 1
BN81 215.1 129.7 129.7 8.6 10.3 140.6 1 1
BN82 214.9 129.8 129.6 7.5 10.3 141.0 1 1
BN83 215.2 129.7 129.1 9.0 9.7 141.9 1 1
BN84 215.2 130.1 129.9 7.9 10.8 141.3 1 1
BN85 215.4 130.7 130.2 9.0 11.1 141.2 1 1
BN86 215.1 129.9 129.6 8.9 10.2 141.5 1 1
BN87 215.2 129.9 129.7 8.7 9.5 141.6 1 1
BN88 215.0 129.6 129.2 8.4 10.2 142.1 1 1
BN89 214.9 130.3 129.9 7.4 11.2 141.5 1 1
BN90 215.0 129.9 129.7 8.0 10.5 142.0 1 1
BN91 214.7 129.7 129.3 8.6 9.6 141.6 1 1
BN92 215.4 130.0 129.9 8.5 9.7 141.4 1 1
BN93 214.9 129.4 129.5 8.2 9.9 141.5 1 1
BN94 214.5 129.5 129.3 7.4 10.7 141.5 1 1
BN95 214.7 129.6 129.5 8.3 10.0 142.0 1 1
BN96 215.6 129.9 129.9 9.0 9.5 141.7 1 1
BN97 215.0 130.4 130.3 9.1 10.2 141.1 1 1
BN98 214.4 129.7 129.5 8.0 10.3 141.2 1 1
BN99 215.1 130.0 129.8 9.1 10.2 141.5 1 1
BN100 214.7 130.0 129.4 7.8 10.0 141.2 1 1
BN101 214.4 130.1 130.3 9.7 11.7 139.8 0 1
BN102 214.9 130.5 130.2 11.0 11.5 139.5 0 1
BN103 214.9 130.3 130.1 8.7 11.7 140.2 0 1
BN104 215.0 130.4 130.6 9.9 10.9 140.3 0 1
BN105 214.7 130.2 130.3 11.8 10.9 139.7 0 1
BN106 215.0 130.2 130.2 10.6 10.7 139.9 0 1
BN107 215.3 130.3 130.1 9.3 12.1 140.2 0 1
BN108 214.8 130.1 130.4 9.8 11.5 139.9 0 1
BN109 215.0 130.2 129.9 10.0 11.9 139.4 0 1
BN110 215.2 130.6 130.8 10.4 11.2 140.3 0 1
BN111 215.2 130.4 130.3 8.0 11.5 139.2 0 1
BN112 215.1 130.5 130.3 10.6 11.5 140.1 0 1
BN113 215.4 130.7 131.1 9.7 11.8 140.6 0 1
BN114 214.9 130.4 129.9 11.4 11.0 139.9 0 1
BN115 215.1 130.3 130.0 10.6 10.8 139.7 0 1
BN116 215.5 130.4 130.0 8.2 11.2 139.2 0 1
BN117 214.7 130.6 130.1 11.8 10.5 139.8 0 1
BN118 214.7 130.4 130.1 12.1 10.4 139.9 0 1
BN119 214.8 130.5 130.2 11.0 11.0 140.0 0 1
BN120 214.4 130.2 129.9 10.1 12.0 139.2 0 1
BN121 214.8 130.3 130.4 10.1 12.1 139.6 0 1
BN122 215.1 130.6 130.3 12.3 10.2 139.6 0 1
BN123 215.3 130.8 131.1 11.6 10.6 140.2 0 1
BN124 215.1 130.7 130.4 10.5 11.2 139.7 0 1
BN125 214.7 130.5 130.5 9.9 10.3 140.1 0 1
BN126 214.9 130.0 130.3 10.2 11.4 139.6 0 1
BN127 215.0 130.4 130.4 9.4 11.6 140.2 0 1
BN128 215.5 130.7 130.3 10.2 11.8 140.0 0 1
BN129 215.1 130.2 130.2 10.1 11.3 140.3 0 1
BN130 214.5 130.2 130.6 9.8 12.1 139.9 0 1
BN131 214.3 130.2 130.0 10.7 10.5 139.8 0 1
BN132 214.5 130.2 129.8 12.3 11.2 139.2 0 1
BN133 214.9 130.5 130.2 10.6 11.5 139.9 0 1
BN134 214.6 130.2 130.4 10.5 11.8 139.7 0 1
BN135 214.2 130.0 130.2 11.0 11.2 139.5 0 1
BN136 214.8 130.1 130.1 11.9 11.1 139.5 0 1
BN137 214.6 129.8 130.2 10.7 11.1 139.4 0 1
BN138 214.9 130.7 130.3 9.3 11.2 138.3 0 1
BN139 214.6 130.4 130.4 11.3 10.8 139.8 0 1
BN140 214.5 130.5 130.2 11.8 10.2 139.6 0 1
BN141 214.8 130.2 130.3 10.0 11.9 139.3 0 1
BN142 214.7 130.0 129.4 10.2 11.0 139.2 0 1
BN143 214.6 130.2 130.4 11.2 10.7 139.9 0 1
BN144 215.0 130.5 130.4 10.6 11.1 139.9 0 1
BN145 214.5 129.8 129.8 11.4 10.0 139.3 0 1
BN146 214.9 130.6 130.4 11.9 10.5 139.8 0 1
BN147 215.0 130.5 130.4 11.4 10.7 139.9 0 1
BN148 215.3 130.6 130.3 9.3 11.3 138.1 0 1
BN149 214.7 130.2 130.1 10.7 11.0 139.4 0 1
BN150 214.9 129.9 130.0 9.9 12.3 139.4 0 1
BN151 214.9 130.3 129.9 11.9 10.6 139.8 0 1
BN152 214.6 129.9 129.7 11.9 10.1 139.0 0 1
BN153 214.6 129.7 129.3 10.4 11.0 139.3 0 1
BN154 214.5 130.1 130.1 12.1 10.3 139.4 0 1
BN155 214.5 130.3 130.0 11.0 11.5 139.5 0 1
BN156 215.1 130.0 130.3 11.6 10.5 139.7 0 1
BN157 214.2 129.7 129.6 10.3 11.4 139.5 0 1
BN158 214.4 130.1 130.0 11.3 10.7 139.2 0 1
BN159 214.8 130.4 130.6 12.5 10.0 139.3 0 1
BN160 214.6 130.6 130.1 8.1 12.1 137.9 0 1
BN161 215.6 130.1 129.7 7.4 12.2 138.4 0 1
BN162 214.9 130.5 130.1 9.9 10.2 138.1 0 1
BN163 214.6 130.1 130.0 11.5 10.6 139.5 0 1
BN164 214.7 130.1 130.2 11.6 10.9 139.1 0 1
BN165 214.3 130.3 130.0 11.4 10.5 139.8 0 1
BN166 215.1 130.3 130.6 10.3 12.0 139.7 0 1
BN167 216.3 130.7 130.4 10.0 10.1 138.8 0 1
BN168 215.6 130.4 130.1 9.6 11.2 138.6 0 1
BN169 214.8 129.9 129.8 9.6 12.0 139.6 0 1
BN170 214.9 130.0 129.9 11.4 10.9 139.7 0 1
BN171 213.9 130.7 130.5 8.7 11.5 137.8 0 1
BN172 214.2 130.6 130.4 12.0 10.2 139.6 0 1
BN173 214.8 130.5 130.3 11.8 10.5 139.4 0 1
BN174 214.8 129.6 130.0 10.4 11.6 139.2 0 1
BN175 214.8 130.1 130.0 11.4 10.5 139.6 0 1
BN176 214.9 130.4 130.2 11.9 10.7 139.0 0 1
BN177 214.3 130.1 130.1 11.6 10.5 139.7 0 1
BN178 214.5 130.4 130.0 9.9 12.0 139.6 0 1
BN179 214.8 130.5 130.3 10.2 12.1 139.1 0 1
BN180 214.5 130.2 130.4 8.2 11.8 137.8 0 1
BN181 215.0 130.4 130.1 11.4 10.7 139.1 0 1
BN182 214.8 130.6 130.6 8.0 11.4 138.7 0 1
BN183 215.0 130.5 130.1 11.0 11.4 139.3 0 1
BN184 214.6 130.5 130.4 10.1 11.4 139.3 0 1
BN185 214.7 130.2 130.1 10.7 11.1 139.5 0 1
BN186 214.7 130.4 130.0 11.5 10.7 139.4 0 1
BN187 214.5 130.4 130.0 8.0 12.2 138.5 0 1
BN188 214.8 130.0 129.7 11.4 10.6 139.2 0 1
BN189 214.8 129.9 130.2 9.6 11.9 139.4 0 1
BN190 214.6 130.3 130.2 12.7 9.1 139.2 0 1
BN191 215.1 130.2 129.8 10.2 12.0 139.4 0 1
BN192 215.4 130.5 130.6 8.8 11.0 138.6 0 1
BN193 214.7 130.3 130.2 10.8 11.1 139.2 0 1
BN194 215.0 130.5 130.3 9.6 11.0 138.5 0 1
BN195 214.9 130.3 130.5 11.6 10.6 139.8 0 1
BN196 215.0 130.4 130.3 9.9 12.1 139.6 0 1
BN197 215.1 130.3 129.9 10.3 11.5 139.7 0 1
BN198 214.8 130.3 130.4 10.6 11.1 140.0 0 1
BN199 214.7 130.7 130.8 11.2 11.2 139.4 0 1
BN200 214.3 129.9 129.9 10.2 11.5 139.6 0 1
;
ENDDATA
CALC:
@SET( 'OROUTE',1); ! Output to window after this many lines in buffer;
@SET( 'WNLINE',10000); ! (27) Max command window lines (Windows only);
! Get index of dependent variable;
@FOR( vardep( k): depndx = k);
! Count number observations in each class;
nobs0 = @SUM( obs( i) | xdat( i, depndx) #eq# 0: 1);
nobs1 = @SUM( obs( i) | xdat( i, depndx) #eq# 1: 1);
! Compute the means: mu1(k) and mu2(k), within group and overall means;
! First get the sums;
@FOR( varexp( k):
mu0( k) = @SUM( obs( i) | xdat( i, depndx) #eq# 0: xdat( i, k));
mu1( k) = @SUM( obs( i) | xdat( i, depndx) #eq# 1: xdat( i, k));
! and now convert to means, overall and within class;
mutot( k) = (mu0( k) + mu1( k))/( nobs0 + nobs1);
mu0( k) = mu0( k)/ nobs0;
mu1( k) = mu1( k)/ nobs1;
);
@WRITE(' Expvar Mean0 Mean1 MeanTot', @NEWLINE(1));
@FOR( varexp( k):
@WRITE(' ', @FORMAT( varexp( k),'9s'),' ', @FORMAT( mu0(k), '9.4f') ,' ', @FORMAT( mu1(k), '9.4f') ,
' ', @FORMAT( mutot(k), '9.4f') , @NEWLINE(1));
);
! Get ready to compute the covariance matrix for all the data;
! Initialize full matrix, so invertible;
@FOR( var( k1):
@FOR( var( k2):
covar( k1, k2) = 0;
);
covar( k1, k1) = 1;
);
! Now actually compute within class covariance.
The -2 in the denominator compensates for estimating mu0 and mu2;
@FOR( varexp( k1):
@FOR( varexp( k2) | k2 #ge# k1:
covar( k1, k2) =
( @SUM( obs( i) | xdat( i, depndx) #eq# 0: ( xdat( i, k1) - mu0( k1))*( xdat( i, k2) - mu0( k2)))
+ @SUM( obs( i) | xdat( i, depndx) #eq# 1: ( xdat( i, k1) - mu1( k1))*( xdat( i, k2) - mu1( k2))))/( nobs0 + nobs1 -2);
covar( k2, k1) = covar( k1, k2); ! It is symmetric;
);
);
! In case we want to see it;
@WRITE( @NEWLINE(1),' Covariance matrix', @NEWLINE( 1));
@FOR( varexp( k1):
@FOR( varexp( k2):
@WRITE(' ', @FORMAT( covar( k1, k2),'9.3f'));
);
@WRITE( @NEWLINE( 1));
);
@WRITE( @NEWLINE( 1));
! Invert the covariance matrix;
covarINV = @INVERSE( covar);
! In case we want to see the inverse;
@WRITE( @NEWLINE(1),' Inverse of Covariance matrix', @NEWLINE( 1));
@FOR( varexp( k1):
@FOR( varexp( k2):
@WRITE(' ', @FORMAT( covarinv( k1, k2),'9.3f'));
);
@WRITE( @NEWLINE( 1));
);
@WRITE( @NEWLINE( 1));
! The beta's are given (by sigma inverse) * (mu2 - mu1);
! Intuitively, if the correlations = 0, then explanatory variable(EV) k gets
a very positive weight if both variance of EV k is small and
outcome type 1 is associate with a large EV k, and
a very negative weight if both variance of EV k is small and
outcome type 0 is assocated with a very large EV k;
@FOR( varexp( k1):
beta( k1) = @SUM( varexp( k2): covarinv( k1, k2)*( mu1( k2) - mu0( k2)));
);
! Compute -score of the mean point;
beta0 = - @SUM( varexp( k): beta( k)*( mu1(k) + mu0( k))/2);
@WRITE(' ExpVar Beta ', @NEWLINE(1));
@WRITE(' Beta0 ', @FORMAT( beta0, '12.4f'), @NEWLINE( 1));
@FOR( varexp( k):
@WRITE( ' ', @FORMAT( varexp( k),'9s'),' ', @FORMAT( beta( k), '12.4f'), @NEWLINE( 1));
);
! Write score of each observation;
@WRITE( @NEWLINE(1),' Obs Type Score Error', @NEWLINE( 1));
@FOR( obs( i):
SCORE( I) =
BETA0 + @SUM( varexp( J): BETA( J)* xdat(I,J));
@IFC( xdat( I, depndx) #LT# 1:
ERR( I) = @SMAX( 0, SCORE( I));
@ELSE
ERR( i) = @SMAX( 0, -SCORE( i));
);
@WRITE( @FORMAT( i, '4.0f'), ' ', xdat( i, depndx), ' ',
@FORMAT( score( i),'10.2f'), ' ', @FORMAT( ERR( i), '10.2f'),@NEWLINE( 1));
);
! Prepare data for a chart of the dependent variable vs. one other variable;
! Get the index of the explanatory variables x and y axis;
@FOR( varpltx( j): pltxndx = j);
@FOR( varplty( j): pltyndx = j);
! Create two sets, based on the category of the dependent variable;
@FOR( obs( i) | xdat( i, depndx) #EQ# 0:
@INSERT( OBS0, i);
xp0(i) = xdat( i, pltxndx);
yp0(i) = xdat( i, pltyndx);
);
@FOR( obs( i) | xdat( i, depndx) #EQ# 1:
@INSERT( OBS1, i);
xp1(i) = xdat( i, pltxndx);
yp1(i) = xdat( i, pltyndx);
);
! Beware when interpreting scatter plots, a red point may hide a green point;
@CHARTSCATTER( 'Discriminant Analysis in 2-space', var( pltxndx), var( pltyndx),
'Type 0 points', xp0, yp0, 'Type 1 points', xp1, yp1) ;
ENDCALC
|
