[system] / trunk / webwork / system / courseScripts / Regression.pm Repository: Repository Listing bbplugincoursesdistsnplrochestersystemwww

# View of /trunk/webwork/system/courseScripts/Regression.pm

Revision 79 - (download) (as text) (annotate)
Mon Jul 2 19:02:06 2001 UTC (11 years, 10 months ago) by maria
File size: 12849 byte(s)
```initial report
```

```    1 #!/usr/bin/perl -w
2
3 package Regression;
4
5 \$VERSION = 0.1;
6
7 use strict;
8
9 ################################################################
10 use constant TINY => 1e-8;
11 use constant DEBUGGING => 0;
12 ################################################################
13 =pod
14
16
17   Regression.pm -     weighted linear regression package (line+plane fitting)
18
20
21 Regression.pm is a multivariate linear regression package.  That is,
22 it estimates the c coefficients for a line-fit of the type
23
24 y= b(0)*x(0) + b(1)*x1 + b(2)*x2 + ... + b(k)*xk
25
26 given a data set of N observations, each with k independent x variables and
27 one y variable.  Naturally, N must be greater than k---and preferably
28 considerably greater.  Any reasonable undergraduate statistics book will
29 explain what a regression is.  Most of the time, the user will provide a
30 constant ('1') as x(0) for each observation in order to allow the regression
31 package to fit an intercept.
32
34
35 If the sample data for (x1, x2, y) includes (\$x1[\$i], \$x2[\$i], \$y[\$i]) for 0<=\$i<=5, type
36
37 \$reg = Regression->new( 3, "y", [ "const", "x1", "x2" ] );
38
39 for(\$i=0; \$i<6; \$i++){
40   \$reg->include( \$y[\$i], [ 1.0, \$x1[\$i], \$x2[\$i] ] );
41 }
42
43 @coeff= \$reg->theta();
44
45 \$b0 = \$coeff[0][0];
46 \$b1 = \$coeff[0][1];
47 \$b2 = \$coeff[0][2];
48
50
51 =head2 Original Algorithm (ALGOL-60):
52
53   W.  M.  Gentleman, University of Waterloo, "Basic Description
54   For Large, Sparse Or Weighted Linear Least Squares Problems
55   (Algorithm AS 75)," Applied Statistics (1974) Vol 23; No. 3
56
58
59 R=Rbar is an upperright triangular matrix, kept in normalized
60 form with implicit 1's on the diagonal.  D is a diagonal scaling
61 matrix.  These correspond to "standard Regression usage" as
62
63                 X' X  = R' D R
64
65 A backsubsitution routine (in thetacov) allows to invert the R
66 matrix (the inverse is upper-right triangular, too!). Call this
67 matrix H, that is H=R^(-1).
68
69     (X' X)^(-1) = [(R' D^(1/2)') (D^(1/2) R)]^(-1)
70     = [ R^-1 D^(-1/2) ] [ R^-1 D^(-1/2) ]'
71
72
73
75
76 This algorithm is the statistical "standard." Insertion of a new
77 observation can be done one obs at any time (WITH A WEIGHT!),
78 and still only takes a low quadratic time.  The storage space
79 requirement is of quadratic order (in the indep variables). A
80 practically infinite number of observations can easily be
81 processed!
82
84
85 Naturally, Gentleman invented this algorithm.  Adaptation by ivo
86 welch. Alan Miller (alan@dmsmelb.mel.dms.CSIRO.AU) pointed out
87 nicer ways to compute the R^2.
88
90
91 =cut
92 ################################################################
93
94
95 #### let's start with handling of missing data ("nan" or "NaN")
96
97 my \$nan= "NaN";
98 sub isNaN {
99   if (\$_[0] !~ /[0-9nan]/) { die "definitely not a number in NaN: '\$_[0]'"; }
100   return (\$_[0]=~ /NaN/i) || (\$_[0] != \$_[0]);
101 }
102
103
104 ################################################################
105
106 =pod
107
109
110 receives the number of variables on each observations (i.e., an integer) and
111 returns the blessed data structure.  Also takes an optional name for this
112 regression to remember, as well as a reference to a k*1 array of names for
113 the X coefficients.
114
115 =cut
116
117 ################################################################
118 sub new {
119   my \$classname= shift(@_);
120   my \$K= shift(@_); # the number of variables
121   my \$regname= shift(@_) || "with no name";
122
123   if (!defined(\$K)) { die "Regression->new needs at least one argument for the number of variables"; }
124   if (\$K<=1) { die "Cannot run a regression without at least two variables."; }
125
126   sub zerovec {
127     my @rv;
128     for (my \$i=0; \$i<=\$_[0]; ++\$i) { \$rv[\$i]=0; }
129     return \@rv;
130   }
131
132   bless {
133    k => \$K,
134    regname => \$regname,
135    xnames => shift(@_),
136
137    # constantly updated
138    n => 0,
139    sse => 0,
140    syy => 0,
141    sy => 0,
142    wghtn => 0,
143    d => zerovec(\$K),
144    thetabar => zerovec(\$K),
145    rbarsize => (\$K+1)*\$K/2+1,
146    rbar => zerovec((\$K+1)*\$K/2+1),
147
148    # other constants
149    neverabort => 0,
150
151    # computed on demand
152    theta => undef,
153    sigmasq => undef,
154    rsq => undef,
155    adjrsq => undef
156   }, \$classname;
157 }
158
159 ################################################################
160
161 =pod
162
164
165 is used for debugging.
166
167 =cut
168
169 ################################################################
170 sub dump {
171   my \$this= \$_[0];
172   print "****************************************************************\n";
173   print "Regression '\$this->{regname}'\n";
174   print "****************************************************************\n";
175   sub print1val {
176     no strict;
177     print "\$_[1](\$_[2])=\t". ((defined(\$_[0]->{ \$_[2] }) ? \$_[0]->{ \$_[2] } : "intentionally undef"));
178
179     my \$ref=\$_[0]->{ \$_[2] };
180
181     if (ref(\$ref) eq 'ARRAY') {
182       my \$arrayref= \$ref;
183       print " \$#\$arrayref+1 elements:\n";
184       if (\$#\$arrayref>30) {
185   print "\t";
186   for(my \$i=0; \$i<\$#\$arrayref+1; ++\$i) { print "\$i='\$arrayref->[\$i]';"; }
187   print "\n";
188       }
189       else {
190   for(my \$i=0; \$i<\$#\$arrayref+1; ++\$i) { print "\t\$i=\t'\$arrayref->[\$i]'\n"; }
191       }
192     }
193     elsif (ref(\$ref) eq 'HASH') {
194       my \$hashref= \$ref;
195       print " ".scalar(keys(%\$hashref))." elements\n";
196       while (my (\$key, \$val) = each(%\$hashref)) {
197   print "\t'\$key'=>'\$val';\n";
198       }
199     }
200     else {
201       print " [was scalar]\n"; }
202   }
203
204   while (my (\$key, \$val) = each(%\$this)) {
205     \$this->print1val(\$key, \$key);
206   }
207   print "****************************************************************\n";
208 }
209
210 ################################################################
211 =pod
212
214
215 prints the estimated coefficients, and R^2 and N.
216
217 =cut
218 ################################################################
219 sub print {
220   my \$this= \$_[0];
221   print "****************************************************************\n";
222   print "Regression '\$this->{regname}'\n";
223   print "****************************************************************\n";
224
225   my \$theta= \$this->theta();
226
227   for (my \$i=0; \$i< \$this->k(); ++\$i) {
228     print "Theta[\$i".(defined(\$this->{xnames}->[\$i]) ? "='\$this->{xnames}->[\$i]'":"")."]= ".sprintf("%12.4f", \$theta->[\$i])."\n";
229   }
230   print "R^2= ".sprintf("%.3f", \$this->rsq()).", N= ".\$this->n()."\n";
231   print "****************************************************************\n";
232 }
233
234
235 ################################################################
236 =pod
237
239
240 receives one new observation.  Call is
241
242   \$blessedregr->include( \$yvariable, [ \$x1, \$x2, \$x3 ... \$xk ], 1.0 );
243
244 where 1.0 is an (optional) weight.  Note that inclusion with a
245 weight of -1 can be used to delete an observation.
246
247 The function returns the number of observations so far included.
248
249 =cut
250 ################################################################
251 sub include {
252   my \$this = shift();
253   my \$yelement= shift();
254   my \$xrow= shift();
255   my \$weight= shift() || 1.0;
256
257   # omit observations with missing observations;
258   if (!defined(\$yelement)) { die "Internal Error: yelement is undef"; }
259   if (isNaN(\$yelement)) { return \$this->{n}; }
260
261   my @xcopy;
262   for (my \$i=1; \$i<=\$this->{k}; ++\$i) {
263     if (!defined(\$xrow->[\$i-1])) { die "Internal Error: xrow [ \$i-1 ] is undef"; }
264     if (isNaN(\$xrow->[\$i-1])) { return \$this->{n}; }
265     \$xcopy[\$i]= \$xrow->[\$i-1];
266   }
267
268   \$this->{syy}+= (\$weight*(\$yelement*\$yelement));
269   \$this->{sy}+= (\$weight*(\$yelement));
270   if (\$weight>=0.0) { ++\$this->{n}; } else { --\$this->{n}; }
271
272   \$this->{wghtn}+= \$weight;
273
274   for (my \$i=1; \$i<=\$this->{k};++\$i) {
275     if (\$weight==0.0) { return \$this->{n}; }
276     if (abs(\$xcopy[\$i])>(TINY)) {
277       my \$xi=\$xcopy[\$i];
278
279       my \$di=\$this->{d}->[\$i];
280       my \$dprimei=\$di+\$weight*(\$xi*\$xi);
281       my \$cbar= \$di/\$dprimei;
282       my \$sbar= \$weight*\$xi/\$dprimei;
283       \$weight*=(\$cbar);
284       \$this->{d}->[\$i]=\$dprimei;
285       my \$nextr=int( ((\$i-1)*( (2.0*\$this->{k}-\$i))/2.0+1) );
286       if (!(\$nextr<=\$this->{rbarsize}) ) { die "Internal Error 2"; }
287       my \$xk;
288       for (my \$kc=\$i+1;\$kc<=\$this->{k};++\$kc) {
289   \$xk=\$xcopy[\$kc]; \$xcopy[\$kc]=\$xk-\$xi*\$this->{rbar}->[\$nextr];
290   \$this->{rbar}->[\$nextr]= \$cbar * \$this->{rbar}->[\$nextr]+\$sbar*\$xk;
291   ++\$nextr;
292       }
293       \$xk=\$yelement; \$yelement-= \$xi*\$this->{thetabar}->[\$i];
294       \$this->{thetabar}->[\$i]= \$cbar*\$this->{thetabar}->[\$i]+\$sbar*\$xk;
295     }
296   }
297   \$this->{sse}+=\$weight*(\$yelement*\$yelement);
298
299   # indicate that Theta is garbage now
300   \$this->{theta}= undef;
301   \$this->{sigmasq}= undef; \$this->{rsq}= undef; \$this->{adjrsq}= undef;
302
303   return \$this->{n};
304 }
305
306
307
308 ################################################################
309
310 =pod
311
313
314 estimates and returns the vector of coefficients.
315
316 =cut
317 ################################################################
318
319 sub theta {
320   my \$this= shift();
321
322   if (defined(\$this->{theta})) { return \$this->{theta}; }
323
324   if (\$this->{n} < \$this->{k}) { return undef; }
325   for (my \$i=(\$this->{k}); \$i>=1; --\$i) {
326     \$this->{theta}->[\$i]= \$this->{thetabar}->[\$i];
327     my \$nextr= int ((\$i-1)*((2.0*\$this->{k}-\$i))/2.0+1);
328     if (!(\$nextr<=\$this->{rbarsize})) { die "Internal Error 3"; }
329     for (my \$kc=\$i+1;\$kc<=\$this->{k};++\$kc) {
330       \$this->{theta}->[\$i]-=(\$this->{rbar}->[\$nextr]*\$this->{theta}->[\$kc]);
331       ++\$nextr;
332     }
333   }
334
335   my \$ref= \$this->{theta}; shift(@\$ref); # we are counting from 0
336
337   # if in a scalar context, otherwise please return the array directly
338   return \$this->{theta};
339 }
340
341 ################################################################
342 =pod
343
344 =head2 rsq, adjrsq, sigmasq, ybar, sst, k, n
345
346 These functions provide common auxiliary information.  rsq, adjrsq,
347 sigmasq, sst, and ybar have not been checked but are likely correct.
348 The results are stored for later usage, although this is somewhat
349 unnecessary because the computation is so simple anyway.
350
351 =cut
352
353 ################################################################
354
355 sub rsq {
356   my \$this= shift();
357   return \$this->{rsq}= 1.0- \$this->{sse} / \$this->sst();
358 }
359
360 sub adjrsq {
361   my \$this= shift();
362   return \$this->{adjrsq}= 1.0- (1.0- \$this->rsq())*(\$this->{n}-1)/(\$this->{n} - \$this->{k});
363 }
364
365 sub sigmasq {
366   my \$this= shift();
367   return \$this->{sigmasq}= (\$this->{n}<=\$this->{k}) ? "Inf" : (\$this->{sse}/(\$this->{n} - \$this->{k}));
368 }
369
370 sub ybar {
371   my \$this= shift();
372   return \$this->{ybar}= \$this->{sy}/\$this->{wghtn};
373 }
374
375 sub sst {
376   my \$this= shift();
377   return \$this->{sst}= (\$this->{syy} - \$this->{wghtn}*(\$this->ybar())**2);
378 }
379
380 sub k {
381   my \$this= shift();
382   return \$this->{k};
383 }
384 sub n {
385   my \$this= shift();
386   return \$this->{n};
387 }
388
389
390 ################################################################
391 =pod
392
393 =head1 DEBUGGING = SAMPLE USAGE CODE
394
395 The sample code included with this package demonstrates regression usage.
396 To execute it, just set the constant DEBUGGING at the script head to 1, and
397 do
398
399   perl Regression.pm
400
401 The printout should be
402
403   ****************************************************************
404   Regression 'sample regression'
405   ****************************************************************
406   Theta[0='const']=       0.2950
407   Theta[1='someX']=       0.6723
408   Theta[2='someY']=       1.0688
409   R^2= 0.808, N= 4
410   ****************************************************************
411
412 =cut
413 ################################################################
414
415 if (DEBUGGING) {
416   package main;
417
418   my \$reg= Statistics::Regression->new( 3, "sample regression", [ "const", "someX", "someY" ] );
419   \$reg->include( 2.0, [ 1.0, 3.0, -1.0 ] );
420   \$reg->include( 1.0, [ 1.0, 5.0, 2.0 ] );
421   \$reg->include( 20.0, [ 1.0, 31.0, 0.0 ] );
422   \$reg->include( 15.0, [ 1.0, 11.0, 2.0 ] );
423
424 #  \$reg->print();   or: my \$coefs= \$reg->theta(); print @coefs; print \$reg->rsq;
425 # my \$coefs= \$reg->theta(); print \$coeff[0];
426 }
427
428 ################################################################
429 =pod
430
432
433 =over 4
434
435 =item Missing
436
437 This package lacks routines to compute the standard errors of
438 the coefficients.  This requires access to a matrix inversion
439 package, and I do not have one at my disposal.  If you want to
440 add one, please let me know.
441
442 =item Perl Problem
443
444 perl is unaware of IEEE number representations.  This makes it a
445 pain to test whether an observation contains any missing
446 variables (coded as 'NaN' in Regression.pm).
447
448 =item Others
449
450 I am a novice perl programmer, so this is probably ugly code.  However, it
451 does seem to work, and I could not find anything equivalent on cpan.
452
453 =back
454
455 =head1 INSTALLATION and DOCUMENTATION
456
457 Installation consists of moving the file 'Regression.pm' into a subdirectory
458 Statistics of your modules path (e.g., /usr/lib/perl5/site_perl/5.6.0/).
459
460 The documentation was produced from the module:
461
462 pod2html -noindex -title "perl weighted least squares regression package" Regression.pm > Regression.html
463
464 The documentation was slightly modified by Maria Voloshina, University of Rochester.
465