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

Annotation of /FigKernelPackages/gjostat.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3