annotate PGAP-1.2.1/Statistics/Distributions.pm @ 11:e0c913351502 draft default tip

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