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

Annotation of /FigKernelPackages/gjostat.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3