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

Diff of /FigKernelPackages/gjostat.pm

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.5, Mon Mar 24 21:08:35 2014 UTC revision 1.6, Fri Aug 29 22:00:24 2014 UTC
# Line 15  Line 15 
15  $mean = mean( @x )  $mean = mean( @x )
16  ( $mean, $stdev ) = mean_stddev( @x )  ( $mean, $stdev ) = mean_stddev( @x )
17  $cc = correl_coef( \@x, \@y )  $cc = correl_coef( \@x, \@y )
18  $cc = correl_coef_z_val( $cc, $n_samples )  $Z  = correl_coef_z_val( $cc, $n_samples )
19    
20    $median = median( @list )
21  $median = general_median( $fraction, @list )  $median = general_median( $fraction, @list )
22    
23  ( $chi_sqr, $df, $n ) = chi_square( \@expected, \@observed )  ( $chi_sqr, $df, $n ) = chi_square( \@expected, \@observed )
# Line 56  Line 57 
57          Should consider explicit multiplication          Should consider explicit multiplication
58    
59          No checking for integer values is performed          No checking for integer values is performed
60    
61  End_of_Usage  End_of_Usage
62  }  }
63    
# Line 69  Line 71 
71          mean_stddev          mean_stddev
72          correl_coef          correl_coef
73          correl_coef_z_val          correl_coef_z_val
74            median
75          general_median          general_median
76          chi_square          chi_square
77          contingency_chi_sqr          contingency_chi_sqr
# Line 112  Line 115 
115  #-----------------------------------------------------------------------------  #-----------------------------------------------------------------------------
116  sub mean  sub mean
117  {  {
118      @_ || return undef;      my $n   = 0;
119      my $sum = 0;      my $sum = 0;
120      foreach ( @_ ) { $sum += $_ }      foreach ( @_ ) { if ( defined $_ ) { $n++; $sum += $_ } }
121      $sum / @_;  
122        $n ? $sum / $n : undef;
123  }  }
124    
125    
# Line 124  Line 128 
128  #-----------------------------------------------------------------------------  #-----------------------------------------------------------------------------
129  sub mean_stddev  sub mean_stddev
130  {  {
131      my $n = @_;      my ( $n, $sum, $sum2 ) = ( 0, 0, 0 );
132      return $n ? ( shift, undef) : ( undef, undef ) if $n < 2;      foreach ( @_ ) { if ( defined $_ ) { $n++; $sum += $_; $sum2 += $_ * $_ } }
133      my ( $sum, $sum2 ) = ( 0, 0 );      return $n ? ( $sum, undef) : ( undef, undef ) if $n < 2;
134      foreach ( @_ ) { $sum += $_; $sum2 += $_ * $_ }  
135      my $x = $sum / $n;      my $x = $sum / $n;
136      ( $x, sqrt( ( $sum2 - $n * $x * $x ) / ( $n - 1 ) ) );      ( $x, sqrt( ( $sum2 - $sum * $x ) / ( $n - 1 ) ) );
137  }  }
138    
139    
# Line 165  Line 169 
169    
170    
171  #-----------------------------------------------------------------------------  #-----------------------------------------------------------------------------
172  #  $cc = correl_coef_z_val( $cc, $n_samples )  #  $Z = correl_coef_z_val( $cc, $n_samples )
173  #  #
174  #  arctanh(x) = ln( ( 1 + x ) / ( 1 - x ) ) / 2  #  arctanh(x) = ln( ( 1 + x ) / ( 1 - x ) ) / 2
175  #-----------------------------------------------------------------------------  #-----------------------------------------------------------------------------
# Line 179  Line 183 
183    
184    
185  #-----------------------------------------------------------------------------  #-----------------------------------------------------------------------------
186    #  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  #  Value that fraction of n sample values are below  #  Value that fraction of n sample values are below
203  #  #
204  #  $median = general_median( $fraction, @list )  #  $median = general_median( $fraction, @list )
205    #
206    #  Musings on splitting the intervals between sample values
207  #-----------------------------------------------------------------------------  #-----------------------------------------------------------------------------
208  #                                     n  #                                     n
209  #        ----------------------------------------------------------------  #        ----------------------------------------------------------------
# Line 195  Line 217 
217  sub general_median  sub general_median
218  {  {
219      my $fr = shift;      my $fr = shift;
220      defined( $fr ) && ( $fr > 0 )      ( defined( $fr ) && ( $fr > 0 ) && ( $fr < 1 ) )
221                     && ( $fr < 1 )          or confess "gjostat::general_median called with bad fraction: $fr\n";
                    || confess "gjostat::general_median called with bad fraction: $fr\n";  
222    
223      my $n = @_;      my @list = sort { $a <=> $b } grep { defined $_ } @_;
224        my $n = @list;
225      my $nbelow = $n * $fr;      my $nbelow = $n * $fr;
226      ( $n > 1 ) && ( $nbelow - 0.5 >=  0 )      ( $n > 1 ) && ( $nbelow - 0.5 >=  0 )
227                 && ( $nbelow + 0.5 <= $n )                 && ( $nbelow + 0.5 <= $n )
228                 || return undef;                 || return undef;
229    
     my (@list) = sort { $a <=> $b } @_;  
230      my $ibelow = int($nbelow - 0.5);      my $ibelow = int($nbelow - 0.5);
231      my $frac   = $nbelow - 0.5 - $ibelow;      my $frac   = $nbelow - 0.5 - $ibelow;
232    
233      my $median = $list[ $ibelow ];      my $median = $list[ $ibelow ];
234      $median += $frac * ( $list[ $ibelow+1 ] - $median )  if $frac;      $median += $frac * ( $list[ $ibelow+1 ] - $median )  if $frac;
235    
236      return $median;      $median;
237  }  }
238    
239    
# Line 706  Line 727 
727      my ($n, $m, $p) = @_;      my ($n, $m, $p) = @_;
728      if ($p == 0) { return ($m ==  0) ? 1 : 0 }      if ($p == 0) { return ($m ==  0) ? 1 : 0 }
729      if ($p == 1) { return ($m == $n) ? 1 : 0 }      if ($p == 1) { return ($m == $n) ? 1 : 0 }
730      return (($n <= 1020) && ($m*log($p) + ($n-$m)*log(1-$p) > -744))      (($n <= 1020) && ($m*log($p) + ($n-$m)*log(1-$p) > -744))
731          ? binomial_coef_0($n, $m) * $p**$m * (1-$p)**($n-$m)          ? 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));          : exp(ln_binomial_coef_0($n, $m) + $m*log($p) + ($n-$m)*log(1-$p));
733  }  }
# Line 1213  Line 1234 
1234      ( 1 == @_ ) || confess "gjostat::factorial called without 1 arg\n";      ( 1 == @_ ) || confess "gjostat::factorial called without 1 arg\n";
1235      my $n = shift;      my $n = shift;
1236      ( $n >= 0 ) && ( $n <= 170)      ( $n >= 0 ) && ( $n <= 170)
1237                  || confess "gjostat::factorial called with out of range argument\n";                  || confess "gjostat::factorial called with out of range argument: '$n'\n";
1238    
1239      factorial_0( $n );      factorial_0( $n );
1240  }  }

Legend:
Removed from v.1.5  
changed lines
  Added in v.1.6

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3