[Bio] / FigKernelPackages / gjostat.pm Repository:
ViewVC logotype

Annotation of /FigKernelPackages/gjostat.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : golsen 1.1 package gjostat;
2 :    
3 : olson 1.5 # This is a SAS component.
4 :    
5 : golsen 1.2 use strict;
6 : golsen 1.3 use Carp qw( confess );
7 : golsen 1.2
8 :     sub usage
9 :     {
10 :     <<'End_of_Usage';
11 :     gjostat.pm - perl functions for dealing with statistics
12 :    
13 :     $var = rand_normal()
14 :    
15 :     $mean = mean( @x )
16 :     ( $mean, $stdev ) = mean_stddev( @x )
17 :     $cc = correl_coef( \@x, \@y )
18 :     $cc = correl_coef_z_val( $cc, $n_samples )
19 :    
20 :     $median = general_median( $fraction, @list )
21 :    
22 :     ( $chi_sqr, $df, $n ) = chi_square( \@expected, \@observed )
23 :     ( $chi_sqr, $df, $n ) = contingency_chi_sqr( @row_refs )
24 :     ( $chi_sqr, $df, $n ) = contingency_chi_sqr( \@row_refs )
25 :     ( $chi_sqr, $df, $n ) = contingency_chi_sqr_2( \@row1, \@row2 )
26 :     ( $chi_sqr, $df, $n ) = contingency_chi_sqr_2( \@row_refs )
27 :    
28 :     $p_value = chisqr_prob( $chisqr, $df )
29 :     $chisqr = chisqr_critical_value( $p_value, $df )
30 :    
31 :     $coef = binomial_coef( $n, $m )
32 :     $prob = binomial_prob_eq_m( $n, $m, $p )
33 :     $prob = binomial_prob_le_m( $n, $m, $p )
34 :     $prob = binomial_prob_ge_m( $n, $m, $p )
35 :    
36 :     $ln_coef = ln_binomial_coef( $n, $m )
37 :     $ln_prob = ln_binomial_prob_eq_m( $n, $m, $p )
38 :    
39 :     $m = binomial_critical_value_m_ge( $n, $p, $P )
40 :     $m = binomial_critical_value_m_le( $n, $p, $P )
41 :    
42 :     $prob = std_normal_ge_z( $z )
43 :     $prob = std_normal_le_z( $z )
44 :    
45 :     $z = std_normal_critical_value_z_ge( $P )
46 :     $z = std_normal_critical_value_z_le( $P )
47 :    
48 :     $prob = poisson_prob_eq_n( $n, $mu )
49 :     $prob = poisson_prob_le_n( $n, $mu )
50 :     $prob = poisson_prob_ge_n( $n, $mu )
51 :    
52 :     $fact = factorial( $n )
53 :     $ln_fact = ln_factorial( $n )
54 :    
55 :     Notes: Exponentiation by int uses **.
56 :     Should consider explicit multiplication
57 :    
58 :     No checking for integer values is performed
59 :     End_of_Usage
60 :     }
61 : golsen 1.1
62 :    
63 :     require Exporter;
64 :    
65 :     our @ISA = qw(Exporter);
66 :     our @EXPORT = qw(
67 :     rand_normal
68 :     mean
69 :     mean_stddev
70 :     correl_coef
71 :     correl_coef_z_val
72 :     general_median
73 :     chi_square
74 :     contingency_chi_sqr
75 :     contingency_chi_sqr_2
76 :     chisqr_prob
77 :     chisqr_critical_value
78 :     binomial_coef
79 :     binomial_prob_eq_m
80 :     binomial_prob_le_m
81 :     binomial_prob_ge_m
82 :     binomial_critical_value_m_ge
83 :     binomial_critical_value_m_le
84 :     ln_binomial_coef
85 :     ln_binomial_prob_eq_m
86 :     std_normal_le_z
87 :     std_normal_ge_z
88 :     std_normal_critical_value_z_ge
89 :     std_normal_critical_value_z_le
90 :     poisson_prob_eq_n
91 :     poisson_prob_le_n
92 :     poisson_prob_ge_n
93 :     factorial
94 :     ln_factorial
95 :     );
96 :    
97 :    
98 :     #-----------------------------------------------------------------------------
99 :     # $var = rand_normal()
100 :     #-----------------------------------------------------------------------------
101 :    
102 :     sub rand_normal
103 :     {
104 :     my $sum = -6;
105 :     for (my $i = 1; $i <= 12; $i++) { $sum += rand }
106 :     $sum;
107 :     }
108 :    
109 :    
110 :     #-----------------------------------------------------------------------------
111 :     # $mean = mean( @x )
112 :     #-----------------------------------------------------------------------------
113 :     sub mean
114 :     {
115 :     @_ || return undef;
116 :     my $sum = 0;
117 :     foreach ( @_ ) { $sum += $_ }
118 :     $sum / @_;
119 :     }
120 :    
121 :    
122 :     #-----------------------------------------------------------------------------
123 :     # ($mean, $stdev) = mean_stddev( @x )
124 :     #-----------------------------------------------------------------------------
125 :     sub mean_stddev
126 :     {
127 :     my $n = @_;
128 :     return $n ? ( shift, undef) : ( undef, undef ) if $n < 2;
129 :     my ( $sum, $sum2 ) = ( 0, 0 );
130 :     foreach ( @_ ) { $sum += $_; $sum2 += $_ * $_ }
131 :     my $x = $sum / $n;
132 :     ( $x, sqrt( ( $sum2 - $n * $x * $x ) / ( $n - 1 ) ) );
133 :     }
134 :    
135 :    
136 :     #-----------------------------------------------------------------------------
137 :     # $cc = correl_coef( \@x, \@y )
138 :     #-----------------------------------------------------------------------------
139 :     sub correl_coef
140 :     {
141 :     my ($xref, $yref) = @_;
142 :     $xref && ref $xref eq 'ARRAY' && $yref && ref $yref eq 'ARRAY'
143 : golsen 1.3 or confess "gjostat::correl_coef() called with invalid parameter types.\n";
144 :     (@$xref == @$yref) || confess "gjostat::correl_coef() called with lists of different lengths\n";
145 : golsen 1.1 (@$xref > 2) || return undef;
146 :     my (@x) = @$xref;
147 :     my (@y) = @$yref;
148 :     my $n = @x;
149 :    
150 :     my ($xsum, $x2sum, $ysum, $y2sum, $xysum) = (0) x 5;
151 :     my ($i, $xi, $yi);
152 :    
153 :     for ($i = 1; $i <= $n; $i++)
154 :     {
155 :     $xi = shift @x; $xsum += $xi; $x2sum += $xi*$xi;
156 :     $yi = shift @y; $ysum += $yi; $y2sum += $yi*$yi;
157 :     $xysum += $xi*$yi;
158 :     }
159 :    
160 :     my $xsd = sqrt( ($x2sum - ($xsum*$xsum/$n)) / ($n - 1) );
161 :     my $ysd = sqrt( ($y2sum - ($ysum*$ysum/$n)) / ($n - 1) );
162 : overbeek 1.4 if (($xsd == 0) || ($ysd == 0)) { return undef }
163 : golsen 1.1 ( $xysum - $xsum * $ysum / $n ) / ( $xsd * $ysd * ( $n - 1 ) );
164 :     }
165 :    
166 :    
167 :     #-----------------------------------------------------------------------------
168 :     # $cc = correl_coef_z_val( $cc, $n_samples )
169 :     #
170 :     # arctanh(x) = ln( ( 1 + x ) / ( 1 - x ) ) / 2
171 :     #-----------------------------------------------------------------------------
172 :     sub correl_coef_z_val
173 :     {
174 :     my( $cc, $n_samples ) = @_;
175 :     defined( $cc ) && $cc >= -1 && $cc <= 1 or return undef;
176 :     defined( $n_samples ) && $n_samples > 2 or return undef;
177 :     sqrt( $n_samples - 3 ) * 0.5 * log( ( 1.000000001 + $cc ) / ( 1.000000001 - $cc ) );
178 :     }
179 :    
180 :    
181 :     #-----------------------------------------------------------------------------
182 :     # Value that fraction of n sample values are below
183 :     #
184 :     # $median = general_median( $fraction, @list )
185 :     #-----------------------------------------------------------------------------
186 :     # n
187 :     # ----------------------------------------------------------------
188 :     # fr 2 3 4 5 6 7
189 :     #------------------------------------------------------------------------
190 :     # 0.25 * 1 0*1 2 0*1 2 3 0*1 2 3 4 0 * 2 3 4 5 0 1*2 3 4 5 6
191 :     # 0.50 0*1 0 * 2 0 1*2 3 0 1 * 3 4 0 1 2*3 4 5 0 1 2 * 4 5 6
192 :     # 0.75 0 * 0 1*2 0 1 2*3 0 1 2 3*4 0 1 2 3 * 5 0 1 2 3 4*5 6
193 :     #------------------------------------------------------------------------
194 :    
195 :     sub general_median
196 :     {
197 :     my $fr = shift;
198 :     defined( $fr ) && ( $fr > 0 )
199 :     && ( $fr < 1 )
200 : golsen 1.3 || confess "gjostat::general_median called with bad fraction: $fr\n";
201 : golsen 1.1
202 :     my $n = @_;
203 :     my $nbelow = $n * $fr;
204 :     ( $n > 1 ) && ( $nbelow - 0.5 >= 0 )
205 :     && ( $nbelow + 0.5 <= $n )
206 :     || return undef;
207 :    
208 :     my (@list) = sort { $a <=> $b } @_;
209 :     my $ibelow = int($nbelow - 0.5);
210 :     my $frac = $nbelow - 0.5 - $ibelow;
211 :    
212 :     my $median = $list[ $ibelow ];
213 :     $median += $frac * ( $list[ $ibelow+1 ] - $median ) if $frac;
214 :    
215 :     return $median;
216 :     }
217 :    
218 :    
219 :     #-----------------------------------------------------------------------------
220 :     # Find the chi-square value for an expected value (or frequency) list and
221 :     # an observed list (expected values need not be normalized):
222 :     #
223 :     # ( $chi_sqr, $df, $n ) = chi_square( \@expected, \@observed )
224 :     #
225 :     #-----------------------------------------------------------------------------
226 :    
227 :     sub chi_square
228 :     {
229 :     my ( $expect, $obs ) = @_;
230 :     $expect && ref $expect eq 'ARRAY' && $obs && ref $obs eq 'ARRAY'
231 : golsen 1.3 or confess "gjostat::chi_square() called with invalid parameter types.";
232 : golsen 1.1 ( @$expect > 1 ) && ( @$expect == @$obs ) || return ( 0, 0, 0 );
233 :    
234 :     my ( $sum1, $sum2 ) = ( 0, 0 );
235 :     foreach ( @$expect ) { $sum1 += $_ }
236 :     foreach ( @$obs ) { $sum2 += $_ }
237 :     ( $sum1 > 0 && $sum2 > 0 ) || return ( 0, 0, 0 );
238 :     my $scale = $sum2 / $sum1;
239 :    
240 :     my ( $e, $o );
241 :     my ( $chisqr, $df ) = ( 0, -1 );
242 :     my $i = 0;
243 :     foreach $e ( map { $scale * $_ } @$expect )
244 :     {
245 :     $o = $obs->[ $i++ ];
246 :     if ( $e > 0 )
247 :     {
248 :     $o -= $e;
249 :     $chisqr += $o * $o / $e;
250 :     $df++;
251 :     }
252 :     elsif ( $o > 0 )
253 :     {
254 : golsen 1.3 confess "gjostat::chi_sqr called with invalid expected value\n"
255 : golsen 1.1 }
256 :     }
257 :    
258 :     ( $df > 0 ) ? ( $chisqr, $df, $sum2 ) : ( 0, 0, 0 )
259 :     }
260 :    
261 :    
262 :     #-----------------------------------------------------------------------------
263 :     # Find the chi-square value for a contingency table:
264 :     #
265 :     # ( $chi_sqr, $df, $n ) = contingency_chi_sqr( @row_refs )
266 :     # ( $chi_sqr, $df, $n ) = contingency_chi_sqr( \@row_refs )
267 :     #
268 :     #-----------------------------------------------------------------------------
269 :    
270 :     sub contingency_chi_sqr
271 :     {
272 :     if ( ( @_ == 1 ) && ( ref( $_[0] ) eq "ARRAY" )
273 :     && ( ref( $_[0]->[0] ) eq "ARRAY" )
274 :     ) { @_ = @{$_[0]} }
275 :    
276 :     ( @_ > 1 ) || return (0, 0, 0);
277 :     ref( $_[0] ) eq "ARRAY"
278 : golsen 1.3 || confess "gjostat::contingency_chi_sqr: arguements must be ARRAY references\n";
279 : golsen 1.1
280 :     my $ncol = @{ $_[0] };
281 :     my ( @rows, @csum, @rsum );
282 :     my $sum = 0;
283 :     my $n;
284 :     my $row;
285 :    
286 :     foreach $row ( @_ )
287 :     {
288 :     ref( $row ) eq "ARRAY"
289 : golsen 1.3 || confess "gjostat::contingency_chi_sqr: arguements must be ARRAY references\n";
290 : golsen 1.1 ( @$row == $ncol )
291 : golsen 1.3 || confess "gjostat::contingency_chi_sqr: all rows must have same number of items\n";
292 : golsen 1.1
293 :     my $rsum = 0;
294 :     for (my $i = 0; $i < $ncol; $i++)
295 :     {
296 :     $n = $row->[ $i ];
297 :     $rsum += $n;
298 :     $csum[ $i ] += $n;
299 :     }
300 :     if ( $rsum > 0 )
301 :     {
302 :     push @rows, $row;
303 :     push @rsum, $rsum;
304 :     $sum += $rsum;
305 :     }
306 :     }
307 :    
308 :     $ncol = 0;
309 :     foreach ( @csum ) { $ncol++ if $_ > 0 }
310 :     ( @rows > 1 ) && ( $ncol > 1 ) || return (0, 0, 0);
311 :    
312 :     my $chi_sqr = 0;
313 :     my ( $e, $rsum, $c );
314 :     foreach $row ( @rows )
315 :     {
316 :     $rsum = shift @rsum;
317 :    
318 :     for (my $i = 0; $i < $ncol; $i++)
319 :     {
320 :     ( ( $c = $csum[ $i ] ) > 0 ) || next;
321 :     $e = $rsum * $c / $sum;
322 :     $chi_sqr += ( $row->[ $i ] - $e )**2 / $e;
323 :     }
324 :     }
325 :    
326 :     ( $chi_sqr, ($ncol-1) * (@rows-1), $sum )
327 :     }
328 :    
329 :    
330 :     #-----------------------------------------------------------------------------
331 :     # Find the chi-square value for a 2-row contingency table:
332 :     #
333 :     # ( $chi_sqr, $df, $n ) = contingency_chi_sqr_2( \@row1, \@row2 )
334 :     # ( $chi_sqr, $df, $n ) = contingency_chi_sqr_2( \@row_refs )
335 :     #
336 :     #-----------------------------------------------------------------------------
337 :    
338 :     sub contingency_chi_sqr_2
339 :     {
340 :     if ( ( @_ == 1 ) && ( ref( $_[0] ) eq "ARRAY" )
341 :     && ( ref( $_[0]->[0] ) eq "ARRAY" )
342 :     ) { @_ = @{$_[0]} }
343 :    
344 :     ( @_ == 2 ) && ( ref( $_[0] ) eq "ARRAY" ) && ( ref( $_[1] ) eq "ARRAY" )
345 : golsen 1.3 || confess "gjostat::contingency_chi_sqr_2: requires two array references\n";
346 : golsen 1.1
347 :     my $ncol = @{ $_[0] };
348 :     ( $ncol == @{ $_[1] } )
349 : golsen 1.3 || confess "gjostat::contingency_chi_sqr_2: all rows must have same number of items\n";
350 : golsen 1.1
351 :     my ( @csum, @rsum );
352 :     my $sum = 0;
353 :     my $n;
354 :    
355 :     foreach my $row ( @_ )
356 :     {
357 :     my $rsum = 0;
358 :     for (my $i = 0; $i < $ncol; $i++)
359 :     {
360 :     $n = $row->[ $i ];
361 :     $rsum += $n;
362 :     $csum[ $i ] += $n;
363 :     }
364 :    
365 :     ( $rsum > 0 ) || return (0, 0, 0);
366 :     push @rsum, $rsum;
367 :     $sum += $rsum;
368 :     }
369 :    
370 :     $ncol = 0;
371 :     foreach ( @csum ) { $ncol++ if $_ > 0 }
372 :     ( $ncol > 1 ) || return (0, 0, 0);
373 :    
374 :     my $chi_sqr = 0;
375 :     my ( $e, $rsum, $c );
376 :     foreach my $row ( @_ )
377 :     {
378 :     $rsum = shift @rsum;
379 :    
380 :     for (my $i = 0; $i < $ncol; $i++)
381 :     {
382 :     ( ( $c = $csum[ $i ] ) > 0 ) || next;
383 :     $e = $rsum * $c / $sum;
384 :     $chi_sqr += ( $row->[ $i ] - $e )**2 / $e;
385 :     }
386 :     }
387 :    
388 :     ( $chi_sqr, $ncol-1, $sum )
389 :     }
390 :    
391 :    
392 :     #-----------------------------------------------------------------------------
393 :     # Probability of a chi square value greater than or equal to chisqr
394 :     #
395 :     # Based on:
396 :     #
397 :     # Zelen, M. and Severo, N. C. (1965). Probability functions. In
398 :     # Handbook of Mathematical Functions, Abramowitz, M. and Stegun,
399 :     # I. A., eds. (New York: Dover Publications), pp. 925-995.
400 :     #
401 :     # Programmed in perl by Gary Olsen
402 :     #
403 :     # $p_value = chisqr_prob( $chisqr, $df )
404 :     #-----------------------------------------------------------------------------
405 :    
406 :     sub chisqr_prob
407 :     {
408 :     my ($chisqr, $df) = @_;
409 :     defined($chisqr) && defined($df)
410 : golsen 1.3 || confess "gjostat::chisqr_prob: undefined arguement\n";
411 : golsen 1.1
412 : golsen 1.3 ($chisqr >= 0) || confess "gjostat::chisqr_prob: bad chi square value: $chisqr\n";
413 : golsen 1.1 ($df > 0) && (int($df) == $df)
414 : golsen 1.3 || confess "gjostat::chisqr_prob: bad degrees of freedom: $df\n";
415 : golsen 1.1
416 :     if ($chisqr == 0) { return 1 }
417 :     if ($chisqr - $df - 10*sqrt($df) > 49.0) { return 1e-14 }
418 :    
419 :     my $inverseProb = 1;
420 :     my $delta = 1;
421 :     my $denom = $df;
422 :    
423 :     while ( 1e16 * $delta > $inverseProb )
424 :     {
425 :     $denom += 2;
426 :     $delta *= $chisqr / $denom;
427 :     $inverseProb += $delta;
428 :     }
429 :    
430 :     $denom = $df;
431 :     my $i;
432 :     for ($i = $df - 2; $i > 0; $i -= 2) { $denom *= $i }
433 :    
434 :     $inverseProb *= exp(-0.5 * $chisqr) * $chisqr**int(0.5*($df+1)) / $denom;
435 :    
436 :     # pi/2 = 1.57079632679489661923
437 :    
438 :     if ($df % 2 != 0) { $inverseProb /= sqrt(1.57079632679489661923 * $chisqr) }
439 :    
440 :     1 - $inverseProb;
441 :     }
442 :    
443 :    
444 :     #-----------------------------------------------------------------------------
445 :     # Find the chi square value required for a given critical P-value:
446 :     #
447 :     # $chisqr = chisqr_critical_value( $p_value, $df )
448 :     #-----------------------------------------------------------------------------
449 :    
450 :     sub chisqr_critical_value
451 :     {
452 : golsen 1.3 ( 2 == @_ ) || confess "gjostat::chisqr_critical_value called without 2 args\n";
453 : golsen 1.1 my ($prob, $df) = @_;
454 :    
455 :     ( $prob > 1e-13 ) && ( $prob < 1 )
456 : golsen 1.3 || confess "gjostat::chisqr_critical_value: chi square out of range: $prob\n";
457 : golsen 1.1
458 :     ( $df > 0 ) && ( $df <= 200 )
459 : golsen 1.3 || confess "gjostat::chisqr_critical_value: df out of range: $df\n";
460 : golsen 1.1
461 :     my $chisqr = $df + 3 * sqrt($df);
462 :     my $step;
463 :    
464 :     if (chisqr_prob($chisqr, $df) > $prob)
465 :     {
466 :     $chisqr *= 2.0;
467 :     while (chisqr_prob($chisqr, $df) > $prob) { $chisqr *= 2.0 }
468 :     $step = 0.25 * $chisqr;
469 :     }
470 :     else
471 :     {
472 :     $step = 0.5 * $chisqr;
473 :     }
474 :    
475 :     $chisqr -= $step;
476 :     $step *= 0.5;
477 :    
478 :     while ($step > 1e-9 * $chisqr)
479 :     {
480 :     if (chisqr_prob($chisqr, $df) > $prob) { $chisqr += $step }
481 :     else { $chisqr -= $step }
482 :     $step *= 0.5;
483 :     }
484 :    
485 :     $chisqr;
486 :     }
487 :    
488 :    
489 :     # Binomial distribution probabilities: ======================================
490 :    
491 :     # $coef = binomial_coef($n, $m)
492 :    
493 :     sub binomial_coef
494 :     {
495 : golsen 1.3 ( 2 == @_ ) || confess "gjostat::binomial_coef called without 2 args\n";
496 : golsen 1.1 my ($n, $m) = @_;
497 :     ( $n >= 1 ) && ( $m >= 0 )
498 :     && ( $m <= $n )
499 : golsen 1.3 || confess "gjostat::binomial_coef called with invalid arg values\n";
500 : golsen 1.1 ( 2 * $m <= $n ) ? binomial_coef_1( $n, $m )
501 :     : binomial_coef_1( $n, $n-$m );
502 :     }
503 :    
504 :    
505 :     # $prob = binomial_prob_eq_m($n, $m, $p)
506 :    
507 :     sub binomial_prob_eq_m
508 :     {
509 : golsen 1.3 ( 3 == @_ ) || confess "gjostat::binomial_prob_eq_m called without 3 args\n";
510 : golsen 1.1 my ($n, $m, $p) = @_;
511 :     ( $n >= 1 ) && ( $m >= 0 )
512 :     && ( $m <= $n )
513 :     && ( $p >= 0 )
514 :     && ( $p <= 1 )
515 : golsen 1.3 || confess "gjostat::binomial_prob_eq_m called with invalid arg values\n";
516 : golsen 1.1 if ( $p == 0 ) { return ($m == 0) ? 1 : 0 }
517 :     if ( $p == 1 ) { return ($m == $n) ? 1 : 0 }
518 :    
519 :     # If no underflow is predicted, to exact calc
520 :    
521 :     ( ( $n <= 1020 ) && ( $m * log($p) + ($n-$m) * log(1-$p) > -744 ) )
522 :     ? binomial_coef_0($n, $m) * $p**$m * (1-$p)**($n-$m)
523 :     : exp(ln_binomial_coef_0($n, $m) + $m*log($p) + ($n-$m)*log(1-$p));
524 :     }
525 :    
526 :    
527 :     sub binomial_prob_le_m
528 :     {
529 : golsen 1.3 ( 3 == @_ ) || confess "gjostat::binomial_prob_le_m called without 3 args\n";
530 : golsen 1.1 my ($n, $m, $p) = @_;
531 :     ( $n >= 1 ) && ( $m >= 0 )
532 :     && ( $m <= $n )
533 :     && ( $p >= 0 )
534 :     && ( $p <= 1 )
535 : golsen 1.3 || confess "gjostat::binomial_prob_le_m called with invalid arg values\n";
536 : golsen 1.1
537 :     if ( ( $p == 0 ) || ( $m == $n ) ) { return 1 }
538 :     if ( $p == 1 ) { return 0 } # special case of m == n handled above
539 :    
540 :     # Figure out the most accurate direction to come from
541 :    
542 :     my ($pn, $w);
543 :     $pn = $p * $n;
544 :     $w = 2 * sqrt($pn * (1-$p));
545 :     ( ($m < $pn - $w) || ( (2 * $m <= $n) && ($m < $pn + $w) ) )
546 :     ? binomial_prob_le_m_00($n, $m, $p)
547 :     : 1 - binomial_prob_ge_m_00($n, $m+1, $p);
548 :     }
549 :    
550 :    
551 :     sub binomial_prob_ge_m
552 :     {
553 : golsen 1.3 ( 3 == @_ ) || confess "gjostat::binomial_prob_ge_m called without 3 args\n";
554 : golsen 1.1 my ($n, $m, $p) = @_;
555 :     ( $n >= 1 ) && ( $m >= 0 )
556 :     && ( $m <= $n )
557 :     && ( $p >= 0 )
558 :     && ( $p <= 1 )
559 : golsen 1.3 || confess "gjostat::binomial_prob_ge_m called with invalid arg values\n";
560 : golsen 1.1
561 :     if ( ( $p == 1 ) || ( $m == 0 ) ) { return 1 }
562 :     if ( $p == 0 ) { return 0 } # special case of $m == 0 handled above
563 :    
564 :     # Figure out the most accurate direction to come from
565 :    
566 :     my ($pn, $w);
567 :     $pn = $p * $n;
568 :     $w = 2 * sqrt($pn * (1-$p));
569 :     ( ($m > $pn + $w) || ( (2 * $m >= $n) && ($m > $pn - $w) ) )
570 :     ? binomial_prob_ge_m_00($n, $m, $p)
571 :     : 1 - binomial_prob_le_m_00($n, $m-1, $p);
572 :     }
573 :    
574 :    
575 :     sub ln_binomial_coef
576 :     {
577 : golsen 1.3 ( 2 == @_ ) || confess "gjostat::ln_binomial_coef called without 2 args\n";
578 : golsen 1.1 my ($n, $m) = @_;
579 :     ( $n >= 1 ) && ( $m >= 0 )
580 :     && ( $m <= $n )
581 : golsen 1.3 || confess "gjostat::ln_binomial_coef called with invalid arg values\n";
582 : golsen 1.1
583 :     ( 2 * $m <= $n ) ? ln_binomial_coef_1($n, $m)
584 :     : ln_binomial_coef_1($n, $n-$m);
585 :     }
586 :    
587 :    
588 :     sub ln_binomial_prob_eq_m
589 :     {
590 : golsen 1.3 ( 3 == @_ ) || confess "gjostat::ln_binomial_prob_eq_m called without 3 args\n";
591 : golsen 1.1 my ($n, $m, $p) = @_;
592 :     ( $n >= 1 ) && ( $m >= 0 )
593 :     && ( $m <= $n )
594 :     && ( $p >= 0 )
595 :     && ( $p <= 1 )
596 : golsen 1.3 || confess "gjostat::ln_binomial_prob_eq_m called with invalid arg values\n";
597 : golsen 1.1 if ( $p == 0 ) { return ($m == 0) ? 0 : undef }
598 :     if ( $p == 1 ) { return ($m == $n) ? 0 : undef }
599 :    
600 :     ln_binomial_coef_0($n, $m) + $m*log($p) + ($n-$m)*log(1-$p);
601 :     }
602 :    
603 :    
604 :     # value of m such that P(>=m) <= P
605 :    
606 :     sub binomial_critical_value_m_ge
607 :     {
608 : golsen 1.3 ( 3 == @_ ) || confess "gjostat::binomial_critical_value_m_ge called without 3 args\n";
609 : golsen 1.1 my ($n, $p, $P) = @_;
610 :     ( $n >= 1 ) && ( $P > 0 )
611 :     && ( $P <= 1 )
612 :     && ( $p >= 0 )
613 :     && ( $p < 1 )
614 : golsen 1.3 || confess "gjostat::binomial_critical_value_m_ge called with invalid arg values\n";
615 : golsen 1.1 if ( $P == 1 ) { return 0 }
616 :     if ( $p == 0 ) { return 1 }
617 :    
618 :     my ($m, $P_ge_m, $term, $q_over_p);
619 :     $P_ge_m = $term = binomial_prob_eq_m_00($n, $n, $p);
620 :     if ( $P_ge_m > $P ) { return undef }
621 :    
622 :     $m = $n;
623 :     $q_over_p = (1-$p) / $p;
624 :     while ( --$m >= 1 )
625 :     {
626 :     if ($P_ge_m == $P) { return $m }
627 :     $P_ge_m += $term *= (($m+1) / ($n-$m)) * $q_over_p;
628 :     if ($P_ge_m > $P) { return $m+1 }
629 :     }
630 :    
631 :     1;
632 :     }
633 :    
634 :    
635 :     sub binomial_critical_value_m_le
636 :     {
637 : golsen 1.3 ( 3 == @_ ) || confess "gjostat::binomial_critical_value_m_le called without 3 args\n";
638 : golsen 1.1 my ($n, $p, $P) = @_;
639 :     ( $n >= 1 ) && ( $P > 0 )
640 :     && ( $P <= 1 )
641 :     && ( $p > 0 )
642 :     && ( $p <= 1 )
643 : golsen 1.3 || confess "gjostat::binomial_critical_value_m_le called with invalid arg values\n";
644 : golsen 1.1 if ($P == 1) { return $n }
645 :     if ($p == 1) { return $n-1 }
646 :    
647 :     my ($m, $P_le_m, $term, $p_over_q);
648 :     $P_le_m = $term = binomial_prob_eq_m_00($n, 0, $p);
649 :     if ($P_le_m > $P) { return undef }
650 :    
651 :     $m = 0;
652 :     $p_over_q = $p / (1-$p);
653 :     while (++$m < $n)
654 :     {
655 :     if ($P_le_m == $P) { return $m }
656 :     $P_le_m += $term *= (($n-$m+1) / $m) * $p_over_q;
657 :     if ($P_le_m > $P) { return $m-1 }
658 :     }
659 :    
660 :     $n-1;
661 :     }
662 :    
663 :    
664 :     # Binomial probability helper functions: ====================================
665 :    
666 :    
667 :     # Given probability of i-1 out of n, what is probability of i out of n?
668 :    
669 :     sub binomial_next_prob_00
670 :     {
671 :     my ($n, $i, $p_over_q, $prob0) = @_;
672 :    
673 :     $prob0 * ( ($n-$i+1) / $i ) * $p_over_q;
674 :     }
675 :    
676 :    
677 :     # Given probability of i+1 out of n, what is probability of i out of n?
678 :    
679 :     sub binomial_prev_prob_00
680 :     {
681 :     my ($n, $i, $p_over_q, $prob0) = @_;
682 :    
683 :     $prob0 * ( ($i+1) / ($n-$i) ) / $p_over_q;
684 :     }
685 :    
686 :    
687 :     # Some basic functions: ====================================================
688 :    
689 :     sub binomial_coef_0 { # no error checking
690 :     my ($n, $m) = @_;
691 :     ( 2 * $m <= $n ) ? binomial_coef_1($n, $m)
692 :     : binomial_coef_1($n, $n-$m);
693 :     }
694 :    
695 :    
696 :     sub binomial_coef_1
697 :     {
698 :     my ($n, $m) = @_;
699 :     my $c = 1;
700 :     while ($m > 0) { $c *= $n-- / $m-- }
701 :     int($c + 0.5);
702 :     }
703 :    
704 :    
705 :     sub binomial_prob_eq_m_0 { # no error checking
706 :     my ($n, $m, $p) = @_;
707 :     if ($p == 0) { return ($m == 0) ? 1 : 0 }
708 :     if ($p == 1) { return ($m == $n) ? 1 : 0 }
709 :     return (($n <= 1020) && ($m*log($p) + ($n-$m)*log(1-$p) > -744))
710 :     ? binomial_coef_0($n, $m) * $p**$m * (1-$p)**($n-$m)
711 :     : exp(ln_binomial_coef_0($n, $m) + $m*log($p) + ($n-$m)*log(1-$p));
712 :     }
713 :    
714 :    
715 :     # Logarithm-based versions for large n: ====================================
716 :    
717 :     sub ln_binomial_coef_0 { # no error checking
718 :     my ($n, $m) = @_;
719 :     ( 2 * $m <= $n ) ? ln_binomial_coef_1($n, $m)
720 :     : ln_binomial_coef_1($n, $n-$m);
721 :     }
722 :    
723 :    
724 :     sub ln_binomial_coef_1
725 :     {
726 :     my ($n, $m) = @_;
727 :     my $c = 0;
728 :     while ($m > 0) { $c += log($n-- / $m--) }
729 :     $c;
730 :     }
731 :    
732 :    
733 :     sub ln_binomial_prob_eq_m_0 { # no error checking
734 :     my ($n, $m, $p) = @_;
735 :     if ($p == 0) { return ($m == 0) ? 1 : 0 }
736 :     if ($p == 1) { return ($m == $n) ? 1 : 0 }
737 :     ln_binomial_coef_0($n, $m) + $m*log($p) + ($n-$m)*log(1-$p);
738 :     }
739 :    
740 :    
741 :     # version for large n?
742 :     #
743 :     # function binomial_prob_eq_m_x_0(n, m, p) { # no error checking
744 :     # if (p == 0) return (m == 0) ? 1 : 0;
745 :     # if (p == 1) return (m == n) ? 1 : 0;
746 :     # if (m <= 0.5 * n) return binomial_prob_eq_m_x_1(n, m, p);
747 :     # else return binomial_prob_eq_m_x_1(n, n-m, 1-p);
748 :     # }
749 :     #
750 :     # function binomial_prob_eq_m_x_1(n, m, p, c, f) { # m <= n/2
751 :     # if (m == 0) return (1-p)^n;
752 :     # f = p * (1-p)^((n-m)/m);
753 :     # c = 1;
754 :     # while (m > 0) c *= f * n-- / m--;
755 :     # return c;
756 :     # }
757 :    
758 :    
759 :     sub binomial_prob_eq_m_00 { # no checking; no special cases
760 :     my ($n, $m, $p) = @_;
761 :     ( ($n <= 1020) && ($m*log($p) + ($n-$m)*log(1-$p) > -744) )
762 :     ? binomial_coef_0($n, $m) * $p**$m * (1-$p)**($n-$m)
763 :     : exp(ln_binomial_coef_0($n, $m) + $m*log($p) + ($n-$m)*log(1-$p));
764 :     }
765 :    
766 :    
767 :     sub binomial_prob_le_m_00
768 :     {
769 :     my ($n, $m, $p) = @_;
770 :     my ($prob, $term, $q_over_p);
771 :    
772 :     $prob = $term = binomial_prob_eq_m_00($n, $m, $p);
773 :     $q_over_p = (1-$p) / $p;
774 :     while (--$m >= 0 && (1e16 * $term > $prob))
775 :     {
776 :     $prob += $term *= (($m+1) / ($n-$m)) * $q_over_p;
777 :     }
778 :    
779 :     $prob;
780 :    
781 :     # prob = term = (1-p)^n;
782 :     # p_over_q = p / (1-p);
783 :     # for (i = 1; i <= m; i++) prob += term *= ((n-i+1) / i) * p_over_q;
784 :     # return prob;
785 :     }
786 :    
787 :    
788 :     sub binomial_prob_ge_m_00
789 :     {
790 :     my ($n, $m, $p) = @_;
791 :     my ($prob, $term, $p_over_q);
792 :    
793 :     $prob = $term = binomial_prob_eq_m_00($n, $m, $p);
794 :     $p_over_q = $p / (1-$p);
795 :     while (++$m <= $n && (1e16 * $term > $prob))
796 :     {
797 :     $prob += $term *= (($n-$m+1) / $m) * $p_over_q;
798 :     }
799 :    
800 :     $prob;
801 :    
802 :     # prob = term = p^n;
803 :     # q_over_p = (1-p) / p;
804 :     # for (i = n-1; i >= m; i--) prob += term *= ((i+1) / (n-i)) * q_over_p;
805 :     # return prob;
806 :     }
807 :    
808 :    
809 :     # Binomial probability slower cleaner versions for checking results: ========
810 :    
811 :    
812 :     # $P = binomial_prob_le_m_slower( $n, $m, $p )
813 :    
814 :     sub binomial_prob_le_m_slower
815 :     {
816 : golsen 1.3 ( 3 == @_ ) || confess "gjostat::binomial_prob_le_m_slower called without 3 args\n";
817 : golsen 1.1 my ($n, $m, $p) = @_;
818 :     ( $n >= 1 ) && ( $m >= 0 )
819 :     && ( $m <= $n )
820 :     && ( $p >= 0 )
821 :     && ( $p <= 1 )
822 : golsen 1.3 || confess "gjostat::binomial_prob_lbinomial_prob_le_m_slowere_m called with invalid arg values\n";
823 : golsen 1.1
824 :     if ( ( $p == 0 ) || ( $m == $n ) ) { return 1 }
825 :     if ( $p == 1 ) { return 0 } # special case of m == n handled above
826 :    
827 :     my ($prob, $term, $p_over_q, $i);
828 :     $prob = $term = (1-$p)**$n;
829 :     $p_over_q = $p / (1-$p);
830 :     for ($i = 1; $i <= $m; $i++)
831 :     {
832 :     $prob += $term = binomial_next_prob_00($n, $i, $p_over_q, $term);
833 :     }
834 :    
835 :     $prob;
836 :     }
837 :    
838 :    
839 :     # $P = binomial_prob_ge_m_slower( $n, $m, $p )
840 :    
841 :     sub binomial_prob_ge_m_slower
842 :     {
843 : golsen 1.3 ( 3 == @_ ) || confess "gjostat::binomial_prob_ge_m_slower called without 3 args\n";
844 : golsen 1.1 my ($n, $m, $p) = @_;
845 :     ( $n >= 1 ) && ( $m >= 0 )
846 :     && ( $m <= $n )
847 :     && ( $p >= 0 )
848 :     && ( $p <= 1 )
849 : golsen 1.3 || confess "gjostat::binomial_prob_ge_m_slower called with invalid arg values\n";
850 : golsen 1.1
851 :     if ( ( $p == 1 ) || ( $m == 0 ) ) { return 1 }
852 :     if ( $p == 0 ) { return 0 } # special case of $m == 0 handled above
853 :    
854 :     my ($prob, $term, $p_over_q, $i);
855 :     $prob = $term = $p**$n;
856 :     $p_over_q = $p / (1-$p);
857 :     for ($i = $n-1; $i >= $m; $i--)
858 :     {
859 :     $prob += $term = binomial_prev_prob_00($n, $i, $p_over_q, $term);
860 :     }
861 :    
862 :     $prob;
863 :     }
864 :    
865 :    
866 :     # $P = binomial_prob_le_m_slowest( $n, $m, $p )
867 :    
868 :     sub binomial_prob_le_m_slowest
869 :     {
870 : golsen 1.3 ( 3 == @_ ) || confess "gjostat::binomial_prob_le_m_slowest called without 3 args\n";
871 : golsen 1.1 my ($n, $m, $p) = @_;
872 :     ( $n >= 1 ) && ( $m >= 0 )
873 :     && ( $m <= $n )
874 :     && ( $p >= 0 )
875 :     && ( $p <= 1 )
876 : golsen 1.3 || confess "gjostat::binomial_prob_le_m_slowest called with invalid arg values\n";
877 : golsen 1.1
878 :     if ( ( $p == 0 ) || ( $m == $n ) ) { return 1 }
879 :     if ( $p == 1 ) { return 0 } # special case of m == n handled above
880 :    
881 :     my ($prob, $i);
882 :     $prob = (1-$p)**$n;
883 :     for ($i = 1; $i <= $m; $i++) { $prob += binomial_prob_eq_m_00($n, $i, $p) }
884 :    
885 :     $prob;
886 :     }
887 :    
888 :    
889 :     # $P = binomial_prob_ge_m_slowest( $n, $m, $p )
890 :    
891 :     sub binomial_prob_ge_m_slowest
892 :     {
893 : golsen 1.3 ( 3 == @_ ) || confess "gjostat::binomial_prob_ge_m_slowest called without 3 args\n";
894 : golsen 1.1 my ($n, $m, $p) = @_;
895 :     ( $n >= 1 ) && ( $m >= 0 )
896 :     && ( $m <= $n )
897 :     && ( $p >= 0 )
898 :     && ( $p <= 1 )
899 : golsen 1.3 || confess "gjostat::binomial_prob_ge_m_slowest called with invalid arg values\n";
900 : golsen 1.1
901 :     if ( ( $p == 1 ) || ( $m == 0 ) ) { return 1 }
902 :     if ( $p == 0 ) { return 0 } # special case of $m == 0 handled above
903 :    
904 :     my ($prob, $i);
905 :     $prob = $p**$n;
906 :     for ($i = $n-$1; $i >= $m; $i--) { $prob += binomial_prob_eq_m_00($n, $i, $p) }
907 :    
908 :     $prob;
909 :     }
910 :    
911 :    
912 :     # Probability under standard normal curve ==================================
913 :    
914 :     # $P = std_normal_le_z( $z )
915 :    
916 :     sub std_normal_le_z
917 :     {
918 : golsen 1.3 ( 1 == @_ ) || confess "gjostat::std_normal_le_z called without 1 arg\n";
919 : golsen 1.1 my ($z) = @_;
920 :     ($z <= 0) ? std_normal_le_z_1( $z)
921 :     : 1 - std_normal_le_z_1(-$z);
922 :     }
923 :    
924 :    
925 :     # $P = std_normal_ge_z( $z )
926 :    
927 :     sub std_normal_ge_z
928 :     {
929 : golsen 1.3 ( 1 == @_ ) || confess "gjostat::std_normal_ge_z called without 1 arg\n";
930 : golsen 1.1 my ($z) = @_;
931 :     ($z >= 0) ? std_normal_le_z_1(-$z)
932 :     : 1 - std_normal_le_z_1( $z);
933 :     }
934 :    
935 :     # $P std_normal_le_z_1( $z ) # -38.4 <= z <= 0
936 :    
937 :     # pi / 2 = 1.57079632679489661923
938 :     # sqrt(2 / pi) = 0.797884560802865
939 :     # sqrt(2 / pi) / 6 = 0.13298076013381091
940 :    
941 :     sub std_normal_le_z_1
942 :     {
943 :     my $step = 0.003;
944 :     my ($z) = @_;
945 :     if ( $z < -38.4 ) { return 0 }
946 :     my ($y, $y2, $p);
947 :     $p = 0;
948 :     $y = exp(-0.5 * $z * $z);
949 :     while ( ( $z >= -38.4 ) && ( 1e15 * $y > $p ) )
950 :     {
951 :     $z -= $step;
952 :     $y += 4 * exp(-0.5 * $z * $z);
953 :     $z -= $step;
954 :     $p += $y + ( $y2 = exp(-0.5 * $z * $z) );
955 :     $y = $y2;
956 :     }
957 :    
958 :     0.13298076013381091 * $p * $step;
959 :     }
960 :    
961 :    
962 :     # $z = std_normal_critical_value_z_ge( $P )
963 :    
964 :     sub std_normal_critical_value_z_ge
965 :     {
966 : golsen 1.3 ( 1 == @_ ) || confess "gjostat::std_normal_critical_value_z_ge called without 1 arg\n";
967 : golsen 1.1 my ($P) = @_;
968 :     ( $P > 0 ) && ( $P < 1 )
969 : golsen 1.3 || confess "gjostat::std_normal_critical_value_z_ge argument out of range\n";
970 : golsen 1.1 ( $P <= 0.5 ) ? std_normal_critical_value_z_ge_2( $P )
971 :     : -std_normal_critical_value_z_ge_2( 1 - $P )
972 :     }
973 :    
974 :    
975 :     # $z = std_normal_critical_value_z_le( $P )
976 :    
977 :     sub std_normal_critical_value_z_le
978 :     {
979 : golsen 1.3 ( 1 == @_ ) || confess "gjostat::std_normal_critical_value_z_le called without 1 arg\n";
980 : golsen 1.1 my ($P) = @_;
981 :     ( $P > 0 ) && ( $P < 1 )
982 : golsen 1.3 || confess "gjostat::std_normal_critical_value_z_le argument out of range\n";
983 : golsen 1.1
984 :     ( $P <= 0.5 ) ? -std_normal_critical_value_z_ge_2( $P )
985 :     : std_normal_critical_value_z_ge_2( 1 - $P )
986 :     }
987 :    
988 :    
989 :     # $Z = std_normal_critical_value_z_ge_2 ( $P ) # 0 < P <= 0.5
990 :     #
991 :     # Good to 8 or more significant figures in Z
992 :     # (not sure why it does not get even better)
993 :    
994 :     sub std_normal_critical_value_z_ge_2
995 :     {
996 :     my ($P) = @_;
997 :     my $P0 = $P / 0.797884560802865; # = sqrt(2 / pi)
998 :    
999 :     # Get a quick estimate using large integration step from infinity
1000 :    
1001 :     my ($step, $step6, $z, $y0, $y1, $y2, $p, $p0);
1002 :    
1003 :     $step = 0.025;
1004 :     $step6 = $step / 6;
1005 :     $z = 38.6;
1006 :     $p = 0;
1007 :     $y2 = exp(-0.5 * $z * $z);
1008 :    
1009 :     while ( $p < $P0 )
1010 :     {
1011 :     $p0 = $p;
1012 :     $y0 = $y2;
1013 :     $z -= $step;
1014 :     $y1 = exp(-0.5 * $z * $z);
1015 :     $z -= $step;
1016 :     $y2 = exp(-0.5 * $z * $z);
1017 :     $p += ( $y0 + 4 * $y1 + $y2 ) * $step6;
1018 :     }
1019 :     $z += 2 * $step;
1020 :     # printf STDERR "; z1 = %11.8f", $z;
1021 :    
1022 :     # Coefs of x**3, x**2 and x**1 in integral of curve from z to z-prime
1023 :    
1024 :     my $c2 = ( 0.5 * ( $y2 + $y0 ) - $y1 ) / 3;
1025 :     my $c1 = $y1 - 0.25 * ( $y2 + 3 * $y0 );
1026 :     my $c0 = $y0;
1027 :    
1028 :     my $dp;
1029 :    
1030 :     my $frac = 1;
1031 :     my $fracstep = $frac;
1032 :     my $step2 = 0.5 * $step;
1033 :    
1034 :     while ( $fracstep > 0.00001 )
1035 :     {
1036 :     $dp = ( ( ( $c2 * $frac ) + $c1 ) * $frac + $c0 ) * $frac * $step2;
1037 :     $fracstep *= 0.5;
1038 :     if ( $p0 + $dp > $P0 ) { $frac -= $fracstep }
1039 :     else { $frac += $fracstep }
1040 :     }
1041 :     my $z0 = $z - $frac * $step;
1042 :     # printf STDERR "; z2 = %11.8f", $z0;
1043 :    
1044 :     # Find out how much we missed P by
1045 :    
1046 :     my $dp0 = ( $P - std_normal_le_z_1(-$z0) ) / 0.797884560802865;
1047 :     if ( $dp0 == 0 ) { return $z0 }
1048 :    
1049 :     # Figure out how far to move to fix the value of the integral
1050 :    
1051 :     $y0 = exp(-0.5 * $z0 * $z0);
1052 :     $step = 2 * $dp0 / $y0; # linear approximation
1053 :     # printf STDERR "; step = %11.8f", $step;
1054 :    
1055 :     my $z1 = $z0 + $step;
1056 :     my $z2 = $z1 + $step;
1057 :     $y1 = exp(-0.5 * $z1 * $z1);
1058 :     $y2 = exp(-0.5 * $z2 * $z2);
1059 :    
1060 :     # Coefs of x**3, x**2 and x**1 in integral of curve from z to z-prime
1061 :    
1062 :     $c2 = ( 0.5 * ( $y2 + $y0 ) - $y1 ) / 3;
1063 :     $c1 = $y1 - 0.25 * ( $y2 + 3 * $y0 );
1064 :     $c0 = $y0;
1065 :    
1066 :     $frac = 1;
1067 :     $fracstep = $frac;
1068 :     $step2 = 0.5 * abs($step);
1069 :     $dp0 = abs($dp0);
1070 :    
1071 :     while ( $fracstep > 0.00001 )
1072 :     {
1073 :     $dp = ( ( ( $c2 * $frac ) + $c1 ) * $frac + $c0) * $frac * $step2;
1074 :     $fracstep *= 0.5;
1075 :     if ( $dp > $dp0 ) { $frac -= $fracstep }
1076 :     else { $frac += $fracstep }
1077 :     }
1078 :    
1079 :     $z0 - $frac * $step;
1080 :     }
1081 :    
1082 :    
1083 :     # $z = std_normal_critical_value_z_ge_1($P) # 0 < P < 0.5
1084 :     #
1085 :     # Earlier, slower version of critical z-value calculation.
1086 :     # Currently not used.
1087 :    
1088 :     sub std_normal_critical_value_z_ge_1
1089 :     {
1090 :     my ($P) = @_;
1091 :     $P /= 0.13298076013381091; # = sqrt(2 / pi) / 6
1092 :    
1093 :     my $step = 0.01;
1094 :     my $z = 38.6;
1095 :     my ($y, $p, $p0);
1096 :     $p = 0;
1097 :     $y = exp(-0.5 * $z * $z) * $step;
1098 :     while ( $p < $P )
1099 :     {
1100 :     $p0 = $p;
1101 :     $z -= $step;
1102 :     $p += $y + 4 * exp(-0.5 * $z * $z) * $step;
1103 :     $z -= $step;
1104 :     $p += ($y = exp(-0.5 * $z * $z) * $step);
1105 :     }
1106 :    
1107 :     $z += 2 * $step;
1108 :     my $y0 = exp(-0.5 * $z * $z);
1109 :     my $stepstep = $step *= 0.5;
1110 :    
1111 :     while ( $stepstep > 1e-10 )
1112 :     {
1113 :     $y = $step * ( $y0
1114 :     + 4 * exp(-0.5 * ($z - $step) * ($z - $step))
1115 :     + exp(-0.5 * ($z - 2*$step) * ($z - 2*$step))
1116 :     );
1117 :     $stepstep *= 0.5;
1118 :     if ( $p0 + $y > $P ) { $step -= $stepstep }
1119 :     else { $step += $stepstep }
1120 :     }
1121 :    
1122 :     $z - 2*$step;
1123 :     }
1124 :    
1125 :    
1126 :    
1127 :     # Poisson distribution probabilities =======================================
1128 :     # none of these include check for integer arguments
1129 :    
1130 :     sub poisson_prob_eq_n
1131 :     {
1132 : golsen 1.3 ( 2 == @_ ) || confess "gjostat::poisson_prob_eq_n called without 2 args\n";
1133 : golsen 1.1 my ($n, $mu) = @_;
1134 :     ( $n >= 0 ) && ( $mu >= 0 )
1135 : golsen 1.3 || confess "gjostat::poisson_prob_eq_n called with invalid arg values\n";
1136 : golsen 1.1
1137 :     # Figure out the most effective approach
1138 :    
1139 :     ( $mu == 0) ? ( ( $n == 0 ) ? 1 : 0 )
1140 :     : poisson_prob_eq_n_0($n, $mu);
1141 :     }
1142 :    
1143 :    
1144 :     sub poisson_prob_le_n
1145 :     {
1146 : golsen 1.3 ( 2 == @_ ) || confess "gjostat::poisson_prob_le_n called without 2 args\n";
1147 : golsen 1.1 my ($n, $mu) = @_;
1148 :     ( $n >= 0 ) && ( $mu >= 0 )
1149 : golsen 1.3 || confess "gjostat::poisson_prob_le_n called with invalid arg values\n";
1150 : golsen 1.1
1151 :     # Figure out the most effective approach
1152 :    
1153 :     ( $mu == 0 ) ? 1
1154 :     : ( $n == 0 ) ? exp(-$mu)
1155 :     : ( $n <= $mu + 5 * sqrt($mu) ) ? poisson_prob_le_n_0($n, $mu)
1156 :     : 1 - poisson_prob_ge_n_0($n+1, $mu);
1157 :     }
1158 :    
1159 :    
1160 :     sub poisson_prob_ge_n
1161 :     {
1162 : golsen 1.3 ( 2 == @_ ) || confess "gjostat::poisson_prob_ge_n called without 2 args\n";
1163 : golsen 1.1 my ($n, $mu) = @_;
1164 :     ( $n >= 0 ) && ( $mu >= 0 )
1165 : golsen 1.3 || confess "gjostat::poisson_prob_ge_n called with invalid arg values\n";
1166 : golsen 1.1
1167 :     # Figure out the most effective approach
1168 :    
1169 :     ( $n == 0 ) ? 1
1170 :     : ( $mu == 0 ) ? 0 # n == 0 handled above
1171 :     : ( $n >= $mu + 5 * sqrt($mu) ) ? poisson_prob_ge_n_0($n, $mu)
1172 :     : 1 - poisson_prob_le_n_0($n-1, $mu);
1173 :     }
1174 :    
1175 :    
1176 :     sub poisson_prob_le_n_0
1177 :     {
1178 :     my ($n, $mu) = @_;
1179 :     my ($p, $term);
1180 :     $p = $term = poisson_prob_eq_n_0($n, $mu);
1181 :     while ( ( --$n >= 0) && ( 1e16 * $term > $p ) )
1182 :     {
1183 :     $p += ( $term *= ($n+1) / $mu );
1184 :     }
1185 :    
1186 :     $p;
1187 :     }
1188 :    
1189 :    
1190 :     sub poisson_prob_ge_n_0
1191 :     {
1192 :     my ($n, $mu) = @_;
1193 :     my ($p, $term);
1194 :     $p = $term = poisson_prob_eq_n_0($n, $mu);
1195 :     while ( 1e16 * $term > $p ) { $p += ( $term *= $mu / (++$n) ) }
1196 :    
1197 :     $p;
1198 :     }
1199 :    
1200 :    
1201 :     sub poisson_prob_eq_n_0
1202 :     {
1203 :     my ($n, $mu) = @_;
1204 :     if ( $n <= 120 ) { return $mu**$n * exp(-$mu) / factorial_0($n) }
1205 :    
1206 :     my $ln_p = $n * log($mu) - $mu - ln_factorial_0($n);
1207 :     ($ln_p > -744) ? exp($ln_p) : 0;
1208 :     }
1209 :    
1210 :    
1211 :     sub factorial
1212 :     {
1213 : golsen 1.3 ( 1 == @_ ) || confess "gjostat::factorial called without 1 arg\n";
1214 : golsen 1.1 my $n = shift;
1215 :     ( $n >= 0 ) && ( $n <= 170)
1216 : golsen 1.3 || confess "gjostat::factorial called with out of range argument\n";
1217 : golsen 1.1
1218 :     factorial_0( $n );
1219 :     }
1220 :    
1221 :    
1222 :     {
1223 :     my @factorial = (1,1);
1224 :     sub factorial_0
1225 :     {
1226 :     my $n = shift;
1227 :     return $factorial[$n] if @factorial > $n;
1228 :    
1229 :     my $i = @factorial;
1230 :     my $f = $factorial[$i-1];
1231 :     $factorial[$n] = undef;
1232 :     while ( $i <= $n ) { $factorial[$i] = $f *= $i++ }
1233 :     $f;
1234 :     }
1235 :     }
1236 :    
1237 :    
1238 :     sub ln_factorial
1239 :     {
1240 : golsen 1.3 ( 1 == @_ ) || confess "gjostat::ln_factorial called without 1 arg\n";
1241 : golsen 1.1 my $n = shift;
1242 :     ( $n >= 0 ) && ( $n <= 1e7)
1243 : golsen 1.3 || confess "gjostat::ln_factorial called with out of range argument\n";
1244 : golsen 1.1
1245 :     ln_factorial_0( $n );
1246 :     }
1247 :    
1248 :    
1249 :     {
1250 :     my @ln_factorial = (0,0);
1251 :     sub ln_factorial_0
1252 :     {
1253 :     my $n = shift;
1254 :     return $ln_factorial[$n] if @ln_factorial > $n;
1255 :    
1256 :     my $i = @ln_factorial;
1257 :     my $f = $ln_factorial[$i-1];
1258 :     $ln_factorial[$n] = undef;
1259 :     while ( $i <= $n ) { $ln_factorial[$i] = $f += log($i++) }
1260 :     $f;
1261 :     }
1262 :     }
1263 :    
1264 :     # Given probability of i-1, what is probability of i?
1265 :    
1266 :     sub poisson_next_prob
1267 :     {
1268 :     my ($i, $mu, $prob0) = @_;
1269 :     ( $i > 0 ) ? $prob0 * $mu / $i : undef;
1270 :     }
1271 :    
1272 :    
1273 :     # Given probability of i+1, what is probability of i?
1274 :    
1275 :     sub poisson_prev_prob
1276 :     {
1277 :     my ($i, $mu, $prob0) = @_;
1278 :     ( $i >= 0 ) ? $prob0 * ($i+1) / $mu : undef;
1279 :     }
1280 :    
1281 :    
1282 :     1;

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3