annotate PGAP-1.2.1/Statistics/Distributions.pm @ 0:83e62a1aeeeb draft

Uploaded

1 package Statistics::Distributions;
2 
3 use strict;
4 use vars qw($VERSION @ISA @EXPORT @EXPORT_OK);
5 use constant PI => 3.1415926536;
6 use constant SIGNIFICANT => 5; # number of significant digits to be returned
7 
8 require Exporter;
9 
10 @ISA = qw(Exporter AutoLoader);
11 # Items to export into callers namespace by default. Note: do not export
12 # names by default without a very good reason. Use EXPORT_OK instead.
13 # Do not simply export all your public functions/methods/constants.
14 @EXPORT_OK = qw(chisqrdistr tdistr fdistr udistr uprob chisqrprob tprob fprob);
15 $VERSION = '1.02';
16 
17 # Preloaded methods go here.
18    
19 sub chisqrdistr { # Percentage points  X^2(x^2,n)
20 	my ($n, $p) = @_;
21 	if ($n <= 0 || abs($n) - abs(int($n)) != 0) {
22 		die "Invalid n: $n\n"; # degree of freedom
23 	}
24 	if ($p <= 0 || $p > 1) {
25 		die "Invalid p: $p\n"; 
26 	}
27 	return precision_string(_subchisqr($n, $p));
28 }
29 
30 sub udistr { # Percentage points   N(0,1^2)
31 	my ($p) = (@_);
32 	if ($p > 1 || $p <= 0) {
33 		die "Invalid p: $p\n";
34 	}
35 	return precision_string(_subu($p));
36 }
37 
38 sub tdistr { # Percentage points   t(x,n)
39 	my ($n, $p) = @_;
40 	if ($n <= 0 || abs($n) - abs(int($n)) != 0) {
41 		die "Invalid n: $n\n";
42 	}
43 	if ($p <= 0 || $p >= 1) {
44 		die "Invalid p: $p\n";
45 	}
46 	return precision_string(_subt($n, $p));
47 }
48 
49 sub fdistr { # Percentage points  F(x,n1,n2)
50 	my ($n, $m, $p) = @_;
51 	if (($n<=0) || ((abs($n)-(abs(int($n))))!=0)) {
52 		die "Invalid n: $n\n"; # first degree of freedom
53 	}
54 	if (($m<=0) || ((abs($m)-(abs(int($m))))!=0)) {
55 		die "Invalid m: $m\n"; # second degree of freedom
56 	}
57 	if (($p<=0) || ($p>1)) {
58 		die "Invalid p: $p\n";
59 	}
60 	return precision_string(_subf($n, $m, $p));
61 }
62 
63 sub uprob { # Upper probability   N(0,1^2)
64 	my ($x) = @_;
65 	return precision_string(_subuprob($x));
66 }
67 
68 sub chisqrprob { # Upper probability   X^2(x^2,n)
69 	my ($n,$x) = @_;
70 	if (($n <= 0) || ((abs($n) - (abs(int($n)))) != 0)) {
71 		die "Invalid n: $n\n"; # degree of freedom
72 	}
73 	return precision_string(_subchisqrprob($n, $x));
74 }
75 
76 sub tprob { # Upper probability   t(x,n)
77 	my ($n, $x) = @_;
78 	if (($n <= 0) || ((abs($n) - abs(int($n))) !=0)) {
79 		die "Invalid n: $n\n"; # degree of freedom
80 	}
81 	return precision_string(_subtprob($n, $x));
82 }
83 
84 sub fprob { # Upper probability   F(x,n1,n2)
85 	my ($n, $m, $x) = @_;
86 	if (($n<=0) || ((abs($n)-(abs(int($n))))!=0)) {
87 		die "Invalid n: $n\n"; # first degree of freedom
88 	}
89 	if (($m<=0) || ((abs($m)-(abs(int($m))))!=0)) {
90 		die "Invalid m: $m\n"; # second degree of freedom
91 	} 
92 	return precision_string(_subfprob($n, $m, $x));
93 }
94 
95 
96 sub _subfprob {
97 	my ($n, $m, $x) = @_;
98 	my $p;
99 
00 	if ($x<=0) {
01 		$p=1;
02 	} elsif ($m % 2 == 0) {
03 		my $z = $m / ($m + $n * $x);
04 		my $a = 1;
05 		for (my $i = $m - 2; $i >= 2; $i -= 2) {
06 			$a = 1 + ($n + $i - 2) / $i * $z * $a;
07 		}
08 		$p = 1 - ((1 - $z) ** ($n / 2) * $a);
09 	} elsif ($n % 2 == 0) {
10 		my $z = $n * $x / ($m + $n * $x);
11 		my $a = 1;
12 		for (my $i = $n - 2; $i >= 2; $i -= 2) {
13 			$a = 1 + ($m + $i - 2) / $i * $z * $a;
14 		}
15 		$p = (1 - $z) ** ($m / 2) * $a;
16 	} else {
17 		my $y = atan2(sqrt($n * $x / $m), 1);
18 		my $z = sin($y) ** 2;
19 		my $a = ($n == 1) ? 0 : 1;
20 		for (my $i = $n - 2; $i >= 3; $i -= 2) {
21 			$a = 1 + ($m + $i - 2) / $i * $z * $a;
22 		} 
23 		my $b = PI;
24 		for (my $i = 2; $i <= $m - 1; $i += 2) {
25 			$b *= ($i - 1) / $i;
26 		}
27 		my $p1 = 2 / $b * sin($y) * cos($y) ** $m * $a;
28 
29 		$z = cos($y) ** 2;
30 		$a = ($m == 1) ? 0 : 1;
31 		for (my $i = $m-2; $i >= 3; $i -= 2) {
32 			$a = 1 + ($i - 1) / $i * $z * $a;
33 		}
34 		$p = max(0, $p1 + 1 - 2 * $y / PI
35 			- 2 / PI * sin($y) * cos($y) * $a);
36 	}
37 	return $p;
38 }
39 
40 
41 sub _subchisqrprob {
42 	my ($n,$x) = @_;
43 	my $p;
44 
45 	if ($x <= 0) {
46 		$p = 1;
47 	} elsif ($n > 100) {
48 		$p = _subuprob((($x / $n) ** (1/3)
49 				- (1 - 2/9/$n)) / sqrt(2/9/$n));
50 	} elsif ($x > 400) {
51 		$p = 0;
52 	} else {   
53 		my ($a, $i, $i1);
54 		if (($n % 2) != 0) {
55 			$p = 2 * _subuprob(sqrt($x));
56 			$a = sqrt(2/PI) * exp(-$x/2) / sqrt($x);
57 			$i1 = 1;
58 		} else {
59 			$p = $a = exp(-$x/2);
60 			$i1 = 2;
61 		}
62 
63 		for ($i = $i1; $i <= ($n-2); $i += 2) {
64 			$a *= $x / $i;
65 			$p += $a;
66 		}
67 	}
68 	return $p;
69 }
70 
71 sub _subu {
72 	my ($p) = @_;
73 	my $y = -log(4 * $p * (1 - $p));
74 	my $x = sqrt(
75 		$y * (1.570796288
76 		  + $y * (.03706987906
77 		  	+ $y * (-.8364353589E-3
78 			  + $y *(-.2250947176E-3
79 			  	+ $y * (.6841218299E-5
80 				  + $y * (0.5824238515E-5
81 					+ $y * (-.104527497E-5
82 					  + $y * (.8360937017E-7
83 						+ $y * (-.3231081277E-8
84 						  + $y * (.3657763036E-10
85 							+ $y *.6936233982E-12)))))))))));
86 	$x = -$x if ($p>.5);
87 	return $x;
88 }
89 
90 sub _subuprob {
91 	my ($x) = @_;
92 	my $p = 0; # if ($absx > 100)
93 	my $absx = abs($x);
94 
95 	if ($absx < 1.9) {
96 		$p = (1 +
97 			$absx * (.049867347
98 			  + $absx * (.0211410061
99 			  	+ $absx * (.0032776263
00 				  + $absx * (.0000380036
01 					+ $absx * (.0000488906
02 					  + $absx * .000005383)))))) ** -16/2;
03 	} elsif ($absx <= 100) {
04 		for (my $i = 18; $i >= 1; $i--) {
05 			$p = $i / ($absx + $p);
06 		}
07 		$p = exp(-.5 * $absx * $absx) 
08 			/ sqrt(2 * PI) / ($absx + $p);
09 	}
10 
11 	$p = 1 - $p if ($x<0);
12 	return $p;
13 }
14 
15    
16 sub _subt {
17 	my ($n, $p) = @_;
18 
19 	if ($p >= 1 || $p <= 0) {
20 		die "Invalid p: $p\n";
21 	}
22 
23 	if ($p == 0.5) {
24 		return 0;
25 	} elsif ($p < 0.5) {
26 		return - _subt($n, 1 - $p);
27 	}
28 
29 	my $u = _subu($p);
30 	my $u2 = $u ** 2;
31 
32 	my $a = ($u2 + 1) / 4;
33 	my $b = ((5 * $u2 + 16) * $u2 + 3) / 96;
34 	my $c = (((3 * $u2 + 19) * $u2 + 17) * $u2 - 15) / 384;
35 	my $d = ((((79 * $u2 + 776) * $u2 + 1482) * $u2 - 1920) * $u2 - 945) 
36 				/ 92160;
37 	my $e = (((((27 * $u2 + 339) * $u2 + 930) * $u2 - 1782) * $u2 - 765) * $u2
38 			+ 17955) / 368640;
39 
40 	my $x = $u * (1 + ($a + ($b + ($c + ($d + $e / $n) / $n) / $n) / $n) / $n);
41 
42 	if ($n <= log10($p) ** 2 + 3) {
43 		my $round;
44 		do { 
45 			my $p1 = _subtprob($n, $x);
46 			my $n1 = $n + 1;
47 			my $delta = ($p1 - $p) 
48 				/ exp(($n1 * log($n1 / ($n + $x * $x)) 
49 					+ log($n/$n1/2/PI) - 1 
50 					+ (1/$n1 - 1/$n) / 6) / 2);
51 			$x += $delta;
52 			$round = sprintf("%.".abs(int(log10(abs $x)-4))."f",$delta);
53 		} while (($x) && ($round != 0));
54 	}
55 	return $x;
56 }
57 
58 sub _subtprob {
59 	my ($n, $x) = @_;
60 
61 	my ($a,$b);
62 	my $w = atan2($x / sqrt($n), 1);
63 	my $z = cos($w) ** 2;
64 	my $y = 1;
65 
66 	for (my $i = $n-2; $i >= 2; $i -= 2) {
67 		$y = 1 + ($i-1) / $i * $z * $y;
68 	} 
69 
70 	if ($n % 2 == 0) {
71 		$a = sin($w)/2;
72 		$b = .5;
73 	} else {
74 		$a = ($n == 1) ? 0 : sin($w)*cos($w)/PI;
75 		$b= .5 + $w/PI;
76 	}
77 	return max(0, 1 - $b - $a * $y);
78 }
79 
80 sub _subf {
81 	my ($n, $m, $p) = @_;
82 	my $x;
83 
84 	if ($p >= 1 || $p <= 0) {
85 		die "Invalid p: $p\n";
86 	}
87 
88 	if ($p == 1) {
89 		$x = 0;
90 	} elsif ($m == 1) {
91 		$x = 1 / (_subt($n, 0.5 - $p / 2) ** 2);
92 	} elsif ($n == 1) {
93 		$x = _subt($m, $p/2) ** 2;
94 	} elsif ($m == 2) {
95 		my $u = _subchisqr($m, 1 - $p);
96 		my $a = $m - 2;
97 		$x = 1 / ($u / $m * (1 +
98 			(($u - $a) / 2 +
99 				(((4 * $u - 11 * $a) * $u + $a * (7 * $m - 10)) / 24 +
00 					(((2 * $u - 10 * $a) * $u + $a * (17 * $m - 26)) * $u
01 						- $a * $a * (9 * $m - 6)
02 					)/48/$n
03 				)/$n
04 			)/$n));
05 	} elsif ($n > $m) {
06 		$x = 1 / _subf2($m, $n, 1 - $p)
07 	} else {
08 		$x = _subf2($n, $m, $p)
09 	}
10 	return $x;
11 }
12 
13 sub _subf2 {
14 	my ($n, $m, $p) = @_;
15 	my $u = _subchisqr($n, $p);
16 	my $n2 = $n - 2;
17 	my $x = $u / $n * 
18 		(1 + 
19 			(($u - $n2) / 2 + 
20 				(((4 * $u - 11 * $n2) * $u + $n2 * (7 * $n - 10)) / 24 + 
21 					(((2 * $u - 10 * $n2) * $u + $n2 * (17 * $n - 26)) * $u 
22 						- $n2 * $n2 * (9 * $n - 6)) / 48 / $m) / $m) / $m);
23 	my $delta;
24 	do {
25 		my $z = exp(
26 			(($n+$m) * log(($n+$m) / ($n * $x + $m)) 
27 				+ ($n - 2) * log($x)
28 				+ log($n * $m / ($n+$m))
29 				- log(4 * PI)
30 				- (1/$n  + 1/$m - 1/($n+$m))/6
31 			)/2);
32 		$delta = (_subfprob($n, $m, $x) - $p) / $z;
33 		$x += $delta;
34 	} while (abs($delta)>3e-4);
35 	return $x;
36 }
37 
38 sub _subchisqr {
39 	my ($n, $p) = @_;
40 	my $x;
41 
42 	if (($p > 1) || ($p <= 0)) {
43 		die "Invalid p: $p\n";
44 	} elsif ($p == 1){
45 		$x = 0;
46 	} elsif ($n == 1) {
47 		$x = _subu($p / 2) ** 2;
48 	} elsif ($n == 2) {
49 		$x = -2 * log($p);
50 	} else {
51 		my $u = _subu($p);
52 		my $u2 = $u * $u;
53 
54 		$x = max(0, $n + sqrt(2 * $n) * $u 
55 			+ 2/3 * ($u2 - 1)
56 			+ $u * ($u2 - 7) / 9 / sqrt(2 * $n)
57 			- 2/405 / $n * ($u2 * (3 *$u2 + 7) - 16));
58 
59 		if ($n <= 100) {
60 			my ($x0, $p1, $z);
61 			do {
62 				$x0 = $x;
63 				if ($x < 0) {
64 					$p1 = 1;
65 				} elsif ($n>100) {
66 					$p1 = _subuprob((($x / $n)**(1/3) - (1 - 2/9/$n))
67 						/ sqrt(2/9/$n));
68 				} elsif ($x>400) {
69 					$p1 = 0;
70 				} else {
71 					my ($i0, $a);
72 					if (($n % 2) != 0) {
73 						$p1 = 2 * _subuprob(sqrt($x));
74 						$a = sqrt(2/PI) * exp(-$x/2) / sqrt($x);
75 						$i0 = 1;
76 					} else {
77 						$p1 = $a = exp(-$x/2);
78 						$i0 = 2;
79 					}
80 
81 					for (my $i = $i0; $i <= $n-2; $i += 2) {
82 						$a *= $x / $i;
83 						$p1 += $a;
84 					}
85 				}
86 				$z = exp((($n-1) * log($x/$n) - log(4*PI*$x) 
87 					+ $n - $x - 1/$n/6) / 2);
88 				$x += ($p1 - $p) / $z;
89 				$x = sprintf("%.5f", $x);
90 			} while (($n < 31) && (abs($x0 - $x) > 1e-4));
91 		}
92 	}
93 	return $x;
94 }
95 
96 sub log10 {
97 	my $n = shift;
98 	return log($n) / log(10);
99 }
00  
01 sub max {
02 	my $max = shift;
03 	my $next;
04 	while (@_) {
05 		$next = shift;
06 		$max = $next if ($next > $max);
07 	}	
08 	return $max;
09 }
10 
11 sub min {
12 	my $min = shift;
13 	my $next;
14 	while (@_) {
15 		$next = shift;
16 		$min = $next if ($next < $min);
17 	}	
18 	return $min;
19 }
20 
21 sub precision {
22 	my ($x) = @_;
23 	return abs int(log10(abs $x) - SIGNIFICANT);
24 }
25 
26 sub precision_string {
27 	my ($x) = @_;
28 	if ($x) {
29 		return sprintf "%." . precision($x) . "f", $x;
30 	} else {
31 		return "0";
32 	}
33 }
34 
35 
36 # Autoload methods go after =cut, and are processed by the autosplit program.
37 
38 1;
39 __END__
40 # Below is the stub of documentation for your module. You better edit it!
41 
42 =head1 NAME
43 
44 Statistics::Distributions - Perl module for calculating critical values and upper probabilities of common statistical distributions
45 
46 =head1 SYNOPSIS
47 
48   use Statistics::Distributions;
49 
50   $chis=Statistics::Distributions::chisqrdistr (2,.05);
51   print "Chi-squared-crit (2 degrees of freedom, 95th percentile "
52        ."= 0.05 level) = $chis\n";
53   
54   $u=Statistics::Distributions::udistr (.05);
55   print "u-crit (95th percentile = 0.05 level) = $u\n";
56   
57   $t=Statistics::Distributions::tdistr (1,.005);
58   print "t-crit (1 degree of freedom, 99.5th percentile = 0.005 level) "
59        ."= $t\n";
60   
61   $f=Statistics::Distributions::fdistr (1,3,.01);
62   print "F-crit (1 degree of freedom in numerator, 3 degrees of freedom "
63        ."in denominator, 99th percentile = 0.01 level) = $f\n";
64   
65   $uprob=Statistics::Distributions::uprob (-0.85);
66   print "upper probability of the u distribution (u = -0.85): Q(u) "
67        ."= 1-G(u) = $uprob\n";
68   
69   $chisprob=Statistics::Distributions::chisqrprob (3,6.25);
70   print "upper probability of the chi-square distribution (3 degrees "
71        ."of freedom, chi-squared = 6.25): Q = 1-G = $chisprob\n";
72   
73   $tprob=Statistics::Distributions::tprob (3,6.251);
74   print "upper probability of the t distribution (3 degrees of "
75        ."freedom, t = 6.251): Q = 1-G = $tprob\n";
76   
77   $fprob=Statistics::Distributions::fprob (3,5,.625);
78   print "upper probability of the F distribution (3 degrees of freedom "
79        ."in numerator, 5 degrees of freedom in denominator, F = 6.25): "
80        ."Q = 1-G = $fprob\n";
81 
82 =head1 DESCRIPTION
83 
84 This Perl module calculates percentage points (5 significant digits) of the u (standard normal) distribution, the student's t distribution, the chi-square distribution and the F distribution. It can also calculate the upper probability (5 significant digits) of the u (standard normal), the chi-square, the t and the F distribution.
85 These critical values are needed to perform statistical tests, like the u test, the t test, the F test and the chi-squared test, and to calculate confidence intervals.
86 
87 If you are interested in more precise algorithms you could look at:
88  StatLib: http://lib.stat.cmu.edu/apstat/ ; 
89  Applied Statistics Algorithms by Griffiths, P. and Hill, I.D., Ellis Horwood: Chichester (1985)
90 
91 =head1 BUGS 
92 
93 This final version 1.02 has been released after more than one year without a bug report on the previous version 0.07.
94 Nevertheless, if you find any bugs or oddities, please do inform the author. 
95 
96 =head1 INSTALLATION 
97 
98 See perlmodinstall for information and options on installing Perl modules. 
99 
00 =head1 AVAILABILITY 
01 
02 The latest version of this module is available from the Distribution Perl Archive Network (CPAN). Please visit http://www.cpan.org/ to find a CPAN site near you or see http://www.cpan.org/authors/id/M/MI/MIKEK/ .
03 
04 =head1 AUTHOR
05 
06 Michael Kospach <[e-mail]>
07 
08 Nice formating, simplification and bug repair by Matthias Trautner Kromann <[e-mail]>
09 
10 =head1 COPYRIGHT 
11 
12 Copyright 2003 Michael Kospach. All rights reserved. 
13 
14 This library is free software; you can redistribute it and/or modify it under the same terms as Perl itself. 
15 
16 =head1 SEE ALSO
17 
18 Statistics::ChiSquare, Statistics::Table::t, Statistics::Table::F, perl(1).
19 
20 =cut
21 
22 
23 
24 
25 
26 
27 
28 
29