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

Annotation of /FigKernelPackages/gjostat.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3