[system] / trunk / pg / lib / Distributions.pm Repository:
ViewVC logotype

Annotation of /trunk/pg/lib/Distributions.pm

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6346 - (view) (download) (as text)

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

aubreyja at gmail dot com
ViewVC Help
Powered by ViewVC 1.0.9