math_extension.rb in distribution-0.6.0

- old
+ new

@@ -1,5 +1,17 @@
+# The next few requires eventually probably need to go in their own gem. They're all functions and constants used by
+# GSL-adapted pure Ruby math functions.
+require "distribution/math_extension/chebyshev_series"
+require "distribution/math_extension/erfc"
+require "distribution/math_extension/exponential_integral"
+require "distribution/math_extension/gammastar"
+require "distribution/math_extension/gsl_utilities"
+require "distribution/math_extension/incomplete_gamma"
+require "distribution/math_extension/incomplete_beta"
+require "distribution/math_extension/log_utilities"
+
+
 if RUBY_VERSION<"1.9"
   require 'mathn'
   def Prime.each(upper,&block) 
     @primes=Prime.new
     @primes.each do |prime|
@@ -28,10 +40,12 @@
     B8  = -1.0 / 30.0
     B10 =  5.0 / 66.0
     B12 = -691.0 / 2730.0
     B14 =  7.0 / 6.0
     B16 = -3617.0 / 510.0
+
+
     # From statistics2
     def loggamma(x)
       v = 1.0
       while (x < N)
         v *= x
@@ -201,46 +215,78 @@
     # Source:
     # * http://mathworld.wolfram.com/BetaFunction.html
     def beta(x,y)
       (gamma(x)*gamma(y)).quo(gamma(x+y))
     end
+
+    # Get pure-Ruby logbeta
     def logbeta(x,y)
-      (loggamma(x)+loggamma(y))-loggamma(x+y)
+      Beta.log_beta(x,y).first
     end
+
+    # Log beta function conforming to style of lgamma (returns sign in second array index)
+    def lbeta(x,y)
+      Beta.log_beta(x,y)
+    end
+
     # I_x(a,b): Regularized incomplete beta function
+    # Fast version. For a exact calculation, based on factorial
+    # use exact_regularized_beta_function
+    def regularized_beta(x,a,b)
+      return 1 if x==1
+      IncompleteBeta.evaluate(a,b,x)
+    end
+    # I_x(a,b): Regularized incomplete beta function
     # TODO: Find a faster version.
     # Source:
     # * http://dlmf.nist.gov/8.17
-    def regularized_beta_function(x,a,b)
+    def exact_regularized_beta(x,a,b)
       return 1 if x==1
-      #incomplete_beta(x,a,b).quo(beta(a,b))
       m=a.to_i
       n=(b+a-1).to_i
       (m..n).inject(0) {|sum,j|
         sum+(binomial_coefficient(n,j)* x**j * (1-x)**(n-j))
       }
 
-    end
-    # B_x(a,b) : Incomplete beta function
-    # Should be replaced by
-    # http://lib.stat.cmu.edu/apstat/63
+     end
+    # 
+    # Incomplete beta function: B(x;a,b)
+    # +a+ and +b+ are parameters and +x+ is 
+    # integration upper limit.
     def incomplete_beta(x,a,b)
-      raise "Doesn't work"
+      IncompleteBeta.evaluate(a,b,x)*beta(a,b)
+      #Math::IncompleteBeta.axpy(1.0, 0.0, a,b,x)
     end
     
     # Rising factorial
     def rising_factorial(x,n)
       factorial(x+n-1).quo(factorial(x-1))
     end
    
     # Ln of gamma
     def loggamma(x)
-      lg=Math.lgamma(x)
-      lg[0]*lg[1]
+      Math.lgamma(x).first
     end
+
+    def incomplete_gamma(a, x = 0, with_error = false)
+      IncompleteGamma.p(a,x, with_error)
+    end
+    alias :gammp :incomplete_gamma
+
+    def gammq(a, x, with_error = false)
+      IncompleteGamma.q(a,x,with_error)
+    end
+
+    def unnormalized_incomplete_gamma(a, x, with_error = false)
+      IncompleteGamma.unnormalized(a, x, with_error)
+    end
+
+    # Not the same as erfc. This is the GSL version, which may have slightly different results.
+    def erfc_e x, with_error = false
+      Erfc.evaluate(x, with_error)
+    end
     
-    
     # Sequences without repetition. n^k'
     # Also called 'failing factorial'
     def permutations(n,k)
       return 1 if k==0
       return n if k==1
@@ -304,16 +350,16 @@
   end
 end
 
 module Math
   include Distribution::MathExtension
-  module_function :factorial, :beta, :loggamma, :binomial_coefficient, :binomial_coefficient_gamma, :regularized_beta_function, :incomplete_beta, :permutations, :rising_factorial , :fast_factorial, :combinations, :logbeta
+  module_function :factorial, :beta, :loggamma, :erfc_e, :unnormalized_incomplete_gamma, :incomplete_gamma, :gammp, :gammq, :binomial_coefficient, :binomial_coefficient_gamma, :exact_regularized_beta, :incomplete_beta, :regularized_beta, :permutations, :rising_factorial , :fast_factorial, :combinations, :logbeta, :lbeta
 end
 
 # Necessary on Ruby 1.9
 module CMath # :nodoc:
   include Distribution::MathExtension
-  module_function :factorial, :beta, :loggamma,  :binomial_coefficient, :binomial_coefficient_gamma, :regularized_beta_function, :incomplete_beta, :permutations, :rising_factorial, :fast_factorial, :combinations, :logbeta
+  module_function :factorial, :beta, :loggamma, :unnormalized_incomplete_gamma, :incomplete_gamma, :gammp, :gammq, :erfc_e, :binomial_coefficient, :binomial_coefficient_gamma,  :incomplete_beta, :exact_regularized_beta, :regularized_beta, :permutations, :rising_factorial, :fast_factorial, :combinations, :logbeta, :lbeta
 end
 
 if RUBY_VERSION<"1.9"
   module Math
     remove_method :loggamma