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

Annotation of /FigKernelPackages/gjostat.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3