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

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

Thu Jun 28 16:59:20 2001 UTC (11 years, 10 months ago) by maria
File size: 11736 byte(s)
```modified documentation
```

```    1 #package Statistics::Distributions;
2 package Distributions;
3
4 # In order to use the functions from this module in your WebWork problems, you have to load PGstatisticsmacros.pl
5
6 # The documentation at the end of this file is slightly modified for WebWork by Maria Voloshina, University of Rochester.
7
8 use strict;
9 use vars qw(\$VERSION @ISA @EXPORT @EXPORT_OK);
10 use constant PI => 3.1415926536;
11 use constant SIGNIFICANT => 6; # number of significant digits to be returned
12
13 require Exporter;
14
16 # Items to export into callers namespace by default. Note: do not export
17 # names by default without a very good reason. Use EXPORT_OK instead.
18 # Do not simply export all your public functions/methods/constants.
19 @EXPORT_OK = qw(chisqrdistr tdistr fdistr udistr uprob chisqrprob tprob fprob);
20 \$VERSION = '0.07';
21
22 # Preloaded methods go here.
23
24 sub chisqrdistr { # Percentage points  X^2(x^2,n)
25   my (\$n, \$p) = @_;
26   if (\$n <= 0 || abs(\$n) - abs(int(\$n)) != 0) {
27     die "Invalid n: \$n\n"; # degree of freedom
28   }
29   if (\$p <= 0 || \$p > 1) {
30     die "Invalid p: \$p\n";
31   }
32   return precision_string(_subchisqr(\$n, \$p));
33 }
34
35 sub udistr { # Percentage points   N(0,1^2)
36   my (\$p) = (@_);
37   if (\$p > 1 || \$p <= 0) {
38     die "Invalid p: \$p\n";
39   }
40   return precision_string(_subu(\$p));
41 }
42
43 sub tdistr { # Percentage points   t(x,n)
44   my (\$n, \$p) = @_;
45   if (\$n <= 0 || abs(\$n) - abs(int(\$n)) != 0) {
46     die "Invalid n: \$n\n";
47   }
48   if (\$p <= 0 || \$p >= 1) {
49     die "Invalid p: \$p\n";
50   }
51   return precision_string(_subt(\$n, \$p));
52 }
53
54 sub fdistr { # Percentage points  F(x,n1,n2)
55   my (\$n, \$m, \$p) = @_;
56   if ((\$n<=0) || ((abs(\$n)-(abs(int(\$n))))!=0)) {
57     die "Invalid n: \$n\n"; # first degree of freedom
58   }
59   if ((\$m<=0) || ((abs(\$m)-(abs(int(\$m))))!=0)) {
60     die "Invalid m: \$m\n"; # second degree of freedom
61   }
62   if ((\$p<=0) || (\$p>1)) {
63     die "Invalid p: \$p\n";
64   }
65   return precision_string(_subf(\$n, \$m, \$p));
66 }
67
68 sub uprob { # Upper probability   N(0,1^2)
69   my (\$x) = @_;
70   return precision_string(_subuprob(\$x));
71 }
72
73 sub chisqrprob { # Upper probability   X^2(x^2,n)
74   my (\$n,\$x) = @_;
75   if ((\$n <= 0) || ((abs(\$n) - (abs(int(\$n)))) != 0)) {
76     die "Invalid n: \$n\n"; # degree of freedom
77   }
78   return precision_string(_subchisqrprob(\$n, \$x));
79 }
80
81 sub tprob { # Upper probability   t(x,n)
82   my (\$n, \$x) = @_;
83   if ((\$n <= 0) || ((abs(\$n) - abs(int(\$n))) !=0)) {
84     die "Invalid n: \$n\n"; # degree of freedom
85   }
86   return precision_string(_subtprob(\$n, \$x));
87 }
88
89 sub fprob { # Upper probability   F(x,n1,n2)
90   my (\$n, \$m, \$x) = @_;
91   if ((\$n<=0) || ((abs(\$n)-(abs(int(\$n))))!=0)) {
92     die "Invalid n: \$n\n"; # first degree of freedom
93   }
94   if ((\$m<=0) || ((abs(\$m)-(abs(int(\$m))))!=0)) {
95     die "Invalid m: \$m\n"; # second degree of freedom
96   }
97   return precision_string(_subfprob(\$n, \$m, \$x));
98 }
99
100
101 sub _subfprob {
102   my (\$n, \$m, \$x) = @_;
103   my \$p;
104
105   if (\$x<=0) {
106     \$p=1;
107   } elsif (\$m % 2 == 0) {
108     my \$z = \$m / (\$m + \$n * \$x);
109     my \$a = 1;
110     for (my \$i = \$m - 2; \$i >= 2; \$i -= 2) {
111       \$a = 1 + (\$n + \$i - 2) / \$i * \$z * \$a;
112     }
113     \$p = 1 - ((1 - \$z) ** (\$n / 2) * \$a);
114   } elsif (\$n % 2 == 0) {
115     my \$z = \$n * \$x / (\$m + \$n * \$x);
116     my \$a = 1;
117     for (my \$i = \$n - 2; \$i >= 2; \$i -= 2) {
118       \$a = 1 + (\$m + \$i - 2) / \$i * \$z * \$a;
119     }
120     \$p = (1 - \$z) ** (\$m / 2) * \$a;
121   } else {
122     my \$y = atan2(sqrt(\$n * \$x / \$m), 1);
123     my \$z = sin(\$y) ** 2;
124     my \$a = (\$n == 1) ? 0 : 1;
125     for (my \$i = \$n - 2; \$i >= 3; \$i -= 2) {
126       \$a = 1 + (\$m + \$i - 2) / \$i * \$z * \$a;
127     }
128     my \$b = PI;
129     for (my \$i = 2; \$i <= \$m - 1; \$i += 2) {
130       \$b *= (\$i - 1) / \$i;
131     }
132     my \$p1 = 2 / \$b * sin(\$y) * cos(\$y) ** \$m * \$a;
133
134     \$z = cos(\$y) ** 2;
135     \$a = (\$m == 1) ? 0 : 1;
136     for (my \$i = \$m-2; \$i >= 3; \$i -= 2) {
137       \$a = 1 + (\$i - 1) / \$i * \$z * \$a;
138     }
139     \$p = max(0, \$p1 + 1 - 2 * \$y / PI
140       - 2 / PI * sin(\$y) * cos(\$y) * \$a);
141   }
142   return \$p;
143 }
144
145
146 sub _subchisqrprob {
147   my (\$n,\$x) = @_;
148   my \$p;
149
150   if (\$x <= 0) {
151     \$p = 1;
152   } elsif (\$n > 100) {
153     \$p = _subuprob(((\$x / \$n) ** (1/3)
154         - (1 - 2/9/\$n)) / sqrt(2/9/\$n));
155   } elsif (\$x > 400) {
156     \$p = 0;
157   } else {
158     my (\$a, \$i, \$i1);
159     if ((\$n % 2) != 0) {
160       \$p = 2 * _subuprob(sqrt(\$x));
161       \$a = sqrt(2/PI) * exp(-\$x/2) / sqrt(\$x);
162       \$i1 = 1;
163     } else {
164       \$p = \$a = exp(-\$x/2);
165       \$i1 = 2;
166     }
167
168     for (\$i = \$i1; \$i <= (\$n-2); \$i += 2) {
169       \$a *= \$x / \$i;
170       \$p += \$a;
171     }
172   }
173   return \$p;
174 }
175
176 sub _subu {
177   my (\$p) = @_;
178   my \$y = -log(4 * \$p * (1 - \$p));
179   my \$x = sqrt(
180     \$y * (1.570796288
181       + \$y * (.03706987906
182         + \$y * (-.8364353589E-3
183         + \$y *(-.2250947176E-3
184           + \$y * (.6841218299E-5
185           + \$y * (0.5824238515E-5
186           + \$y * (-.104527497E-5
187             + \$y * (.8360937017E-7
188             + \$y * (-.3231081277E-8
189               + \$y * (.3657763036E-10
190               + \$y *.6936233982E-12)))))))))));
191   \$x = -\$x if (\$p>.5);
192   return \$x;
193 }
194
195 sub _subuprob {
196   my (\$x) = @_;
197   my \$p = 0; # if (\$absx > 100)
198   my \$absx = abs(\$x);
199
200   if (\$absx < 1.9) {
201     \$p = (1 +
202       \$absx * (.049867347
203         + \$absx * (.0211410061
204           + \$absx * (.0032776263
205           + \$absx * (.0000380036
206           + \$absx * (.0000488906
207             + \$absx * .000005383)))))) ** -16/2;
208   } elsif (\$absx <= 100) {
209     for (my \$i = 18; \$i >= 1; \$i--) {
210       \$p = \$i / (\$absx + \$p);
211     }
212     \$p = exp(-.5 * \$absx * \$absx)
213       / sqrt(2 * PI) / (\$absx + \$p);
214   }
215
216   \$p = 1 - \$p if (\$x<0);
217   return \$p;
218 }
219
220
221 sub _subt {
222   my (\$n, \$p) = @_;
223
224   if (\$p >= 1 || \$p <= 0) {
225     die "Invalid p: \$p\n";
226   }
227
228   if (\$p == 0.5) {
229     return 0;
230   } elsif (\$p < 0.5) {
231     return - _subt(\$n, 1 - \$p);
232   }
233
234   my \$u = _subu(\$p);
235   my \$u2 = \$u ** 2;
236
237   my \$a = (\$u2 + 1) / 4;
238   my \$b = ((5 * \$u2 + 16) * \$u2 + 3) / 96;
239   my \$c = (((3 * \$u2 + 19) * \$u2 + 17) * \$u2 - 15) / 384;
240   my \$d = ((((79 * \$u2 + 776) * \$u2 + 1482) * \$u2 - 1920) * \$u2 - 945)
241         / 92160;
242   my \$e = (((((27 * \$u2 + 339) * \$u2 + 930) * \$u2 - 1782) * \$u2 - 765) * \$u2
243       + 17955) / 368640;
244
245   my \$x = \$u * (1 + (\$a + (\$b + (\$c + (\$d + \$e / \$n) / \$n) / \$n) / \$n) / \$n);
246
247   if (\$n <= log10(\$p) ** 2 + 3) {
248     my \$round;
249     do {
250       my \$p1 = _subtprob(\$n, \$x);
251       my \$n1 = \$n + 1;
252       my \$delta = (\$p1 - \$p)
253         / exp((\$n1 * log(\$n1 / (\$n + \$x * \$x))
254           + log(\$n/\$n1/2/PI) - 1
255           + (1/\$n1 - 1/\$n) / 6) / 2);
256       \$x += \$delta;
257       \$round = sprintf("%.".abs(int(log10(abs \$x)-4))."f",\$delta);
258     } while ((\$x) && (\$round != 0));
259   }
260   return \$x;
261 }
262
263 sub _subtprob {
264   my (\$n, \$x) = @_;
265
266   my (\$a,\$b);
267   my \$w = atan2(\$x / sqrt(\$n), 1);
268   my \$z = cos(\$w) ** 2;
269   my \$y = 1;
270
271   for (my \$i = \$n-2; \$i >= 2; \$i -= 2) {
272     \$y = 1 + (\$i-1) / \$i * \$z * \$y;
273   }
274
275   if (\$n % 2 == 0) {
276     \$a = sin(\$w)/2;
277     \$b = .5;
278   } else {
279     \$a = (\$n == 1) ? 0 : sin(\$w)*cos(\$w)/PI;
280     \$b= .5 + \$w/PI;
281   }
282   return max(0, 1 - \$b - \$a * \$y);
283 }
284
285 sub _subf {
286   my (\$n, \$m, \$p) = @_;
287   my \$x;
288
289   if (\$p >= 1 || \$p <= 0) {
290     die "Invalid p: \$p\n";
291   }
292
293   if (\$p == 1) {
294     \$x = 0;
295   } elsif (\$m == 1) {
296     \$x = 1 / (_subt(\$n, 0.5 - \$p / 2) ** 2);
297   } elsif (\$n == 1) {
298     \$x = _subt(\$m, \$p/2) ** 2;
299   } elsif (\$m == 2) {
300     my \$u = _subchisqr(\$m, 1 - \$p);
301     my \$a = \$m - 2;
302     \$x = 1 / (\$u / \$m * (1 +
303       ((\$u - \$a) / 2 +
304         (((4 * \$u - 11 * \$a) * \$u + \$a * (7 * \$m - 10)) / 24 +
305           (((2 * \$u - 10 * \$a) * \$u + \$a * (17 * \$m - 26)) * \$u
306             - \$a * \$a * (9 * \$m - 6)
307           )/48/\$n
308         )/\$n
309       )/\$n));
310   } elsif (\$n > \$m) {
311     \$x = 1 / _subf2(\$m, \$n, 1 - \$p)
312   } else {
313     \$x = _subf2(\$n, \$m, \$p)
314   }
315   return \$x;
316 }
317
318 sub _subf2 {
319   my (\$n, \$m, \$p) = @_;
320   my \$u = _subchisqr(\$n, \$p);
321   my \$n2 = \$n - 2;
322   my \$x = \$u / \$n *
323     (1 +
324       ((\$u - \$n2) / 2 +
325         (((4 * \$u - 11 * \$n2) * \$u + \$n2 * (7 * \$n - 10)) / 24 +
326           (((2 * \$u - 10 * \$n2) * \$u + \$n2 * (17 * \$n - 26)) * \$u
327             - \$n2 * \$n2 * (9 * \$n - 6)) / 48 / \$m) / \$m) / \$m);
328   my \$delta;
329   do {
330     my \$z = exp(
331       ((\$n+\$m) * log((\$n+\$m) / (\$n * \$x + \$m))
332         + (\$n - 2) * log(\$x)
333         + log(\$n * \$m / (\$n+\$m))
334         - log(4 * PI)
335         - (1/\$n  + 1/\$m - 1/(\$n+\$m))/6
336       )/2);
337     \$delta = (_subfprob(\$n, \$m, \$x) - \$p) / \$z;
338     \$x += \$delta;
339   } while (abs(\$delta)>3e-4);
340   return \$x;
341 }
342
343 sub _subchisqr {
344   my (\$n, \$p) = @_;
345   my \$x;
346
347   if ((\$p > 1) || (\$p <= 0)) {
348     die "Invalid p: \$p\n";
349   } elsif (\$p == 1){
350     \$x = 0;
351   } elsif (\$n == 1) {
352     \$x = _subu(\$p / 2) ** 2;
353   } elsif (\$n == 2) {
354     \$x = -2 * log(\$p);
355   } else {
356     my \$u = _subu(\$p);
357     my \$u2 = \$u * \$u;
358
359     \$x = max(0, \$n + sqrt(2 * \$n) * \$u
360       + 2/3 * (\$u2 - 1)
361       + \$u * (\$u2 - 7) / 9 / sqrt(2 * \$n)
362       - 2/405 / \$n * (\$u2 * (3 *\$u2 + 7) - 16));
363
364     if (\$n <= 100) {
365       my (\$x0, \$p1, \$z);
366       do {
367         \$x0 = \$x;
368         if (\$x < 0) {
369           \$p1 = 1;
370         } elsif (\$n>100) {
371           \$p1 = _subuprob(((\$x / \$n)**(1/3) - (1 - 2/9/\$n))
372             / sqrt(2/9/\$n));
373         } elsif (\$x>400) {
374           \$p1 = 0;
375         } else {
376           my (\$i0, \$a);
377           if ((\$n % 2) != 0) {
378             \$p1 = 2 * _subuprob(sqrt(\$x));
379             \$a = sqrt(2/PI) * exp(-\$x/2) / sqrt(\$x);
380             \$i0 = 1;
381           } else {
382             \$p1 = \$a = exp(-\$x/2);
383             \$i0 = 2;
384           }
385
386           for (my \$i = \$i0; \$i <= \$n-2; \$i += 2) {
387             \$a *= \$x / \$i;
388             \$p1 += \$a;
389           }
390         }
391         \$z = exp(((\$n-1) * log(\$x/\$n) - log(4*PI*\$x)
392           + \$n - \$x - 1/\$n/6) / 2);
393         \$x += (\$p1 - \$p) / \$z;
394         \$x = sprintf("%.5f", \$x);
395       } while ((\$n < 31) && (abs(\$x0 - \$x) > 1e-4));
396     }
397   }
398   return \$x;
399 }
400
401 sub log10 {
402   my \$n = shift;
403   return log(\$n) / log(10);
404 }
405
406 sub max {
407   my \$max = shift;
408   my \$next;
409   while (@_) {
410     \$next = shift;
411     \$max = \$next if (\$next > \$max);
412   }
413   return \$max;
414 }
415
416 sub min {
417   my \$min = shift;
418   my \$next;
419   while (@_) {
420     \$next = shift;
421     \$min = \$next if (\$next < \$min);
422   }
423   return \$min;
424 }
425
426 sub precision {
427   my (\$x) = @_;
428   return abs int(log10(abs \$x) - SIGNIFICANT);
429 }
430
431 sub precision_string {
432   my (\$x) = @_;
433   if (\$x) {
434     return sprintf "%." . precision(\$x) . "f", \$x;
435   } else {
436     return "0";
437   }
438 }
439
440
441 # Autoload methods go after =cut, and are processed by the autosplit program.
442
443 1;
444 __END__
445
446 # Below is the stub of documentation for your module.
447
449
450 Statistics::Distributions - Perl module for calculating probabilities and critical values of common statistical distributions
451
453
454   use Statistics::Distributions;
455
456   \$uprob=uprob(-0.85);
457   print "upper probability of the u distribution: Q(u) = 1-G(u) (u=1.43) = \$uprob";
458
459   \$tprob=tprob(3, 6.251);
460   print "upper probability of the t distribution: Q = 1-G (3 degrees of freedom  , t = 6.251) = \$tprob";
461
462   \$chisprob=chisqrprob(3, 6.25);
463   print "upper probability of the chi-square distribution: Q = 1-G (3 degrees of freedom, chi-squared = 6.25) = \$chisprob";
464
465   \$fprob=fprob(3,5, .625);
466   print "upper probability of the F distribution: Q = 1-G (3 degrees of freedom  in numerator, 5 degrees of freedom in denominator, F \$
467
468   \$u=udistr(.05);
469   print "u-crit (95th percentile = 0.05 level) = \$u";
470
471   \$t=tdistr(1, .005);
472   print "t-crit (1 degree of freedom, 99.5th percentile = 0.005 level) =\$t";
473
474   \$chis=chisqrdistr(2, .05);
475   print "Chi-squared-crit (2 degrees of freedom, 95th percentile = 0.05 level) = \$chis";
476
477   \$f=fdistr(1,3, .01);
478   print "F-crit (1 degree of freedom in numerator, 3 degrees of freedom in denominator, 99th percentile = 0.01 level) = \$f";
479
481
482 This Perl module calulates percentage points (6 significant digits) of the u (standard normal) distribution,
483 the student's t distribution, the chi-square distribution and the F distribution.
484
485 It can also calculate the upper probability (6 significant digits) of the u (standard normal),
486 the chi-square, the t and the F distribution.
487
488 These critical values are needed to perform statistical tests, like the u test, the t test, the chi-squared test,
489 and the F test, and to calculate confidence intervals.
490
491 If you are interested in more precise algorithms you could look at:
492  StatLib: http://lib.stat.cmu.edu/apstat/
493  Applied Statistics Algorithms by Griffiths, P. and Hill, I.D., Ellis Horwood: Chichester (1985)
494
496
497 Michael Kospach, mike.perl@gmx.at
498 Nice formating, simplification and bug repair by Matthias Trautner Kromann, mtk@id.cbs.dk
499
500 =cut
501
502
503
504
505
506
507
508
509
```