Module: Bullshit::Functions

Extended by:
Math
Included in:
ChiSquareDistribution, NormalDistribution, TDistribution
Defined in:
lib/bullshit.rb

Constant Summary collapse

LANCZOS_COEFFICIENTS = 
[
  0.99999999999999709182,
  57.156235665862923517,
  -59.597960355475491248,
  14.136097974741747174,
  -0.49191381609762019978,
  0.33994649984811888699e-4,
  0.46523628927048575665e-4,
  -0.98374475304879564677e-4,
  0.15808870322491248884e-3,
  -0.21026444172410488319e-3,
  0.21743961811521264320e-3,
  -0.16431810653676389022e-3,
  0.84418223983852743293e-4,
  -0.26190838401581408670e-4,
  0.36899182659531622704e-5,
]
HALF_LOG_2_PI = 
0.5 * log(2 * Math::PI)
ROOT2 = 
sqrt(2)
A = 
-8 * (Math::PI - 3) / (3 * Math::PI * (Math::PI - 4))

Class Method Summary collapse

Instance Method Summary collapse

Class Method Details

.beta_regularized(x, a, b, epsilon = 1E-16, max_iterations = 1 << 16) ⇒ Object

Return an approximation value of Euler’s regularized beta function for x, a, and b with an error <= epsilon, but only iterate max_iterations-times.



750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
 
# File 'lib/bullshit.rb', line 750

def beta_regularized(x, a, b, epsilon = 1E-16, max_iterations = 1 << 16)
  x, a, b = x.to_f, a.to_f, b.to_f
  case
  when a.nan? || b.nan? || x.nan? || a <= 0 || b <= 0 || x < 0 || x > 1
    0 / 0.0
  when x > (a + 1) / (a + b + 2)
    1 - beta_regularized(1 - x, b, a, epsilon, max_iterations)
  else
    fraction = ContinuedFraction.for_b do |n, x|
      if n % 2 == 0
        m = n / 2.0
        (m * (b - m) * x) / ((a + (2 * m) - 1) * (a + (2 * m)))
      else
        m = (n - 1) / 2.0
        -((a + m) * (a + b + m) * x) / ((a + 2 * m) * (a + 2 * m + 1))
      end
    end
    exp(a * log(x) + b * log(1.0 - x) - log(a) - log_beta(a, b)) / 
      fraction[x, epsilon, max_iterations]
  end
rescue Errno::ERANGE, Errno::EDOM
  0 / 0.0
end

.gammaP_regularized(x, a, epsilon = 1E-16, max_iterations = 1 << 16) ⇒ Object

Return an approximation of the regularized gammaP function for x and a with an error of <= epsilon, but only iterate max_iterations-times.



777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
 
# File 'lib/bullshit.rb', line 777

def gammaP_regularized(x, a, epsilon = 1E-16, max_iterations = 1 << 16)
  x, a = x.to_f, a.to_f
  case
  when a.nan? || x.nan? || a <= 0 || x < 0
    0 / 0.0
  when x == 0
    0.0
  when 1 <= a && a < x
    1 - gammaQ_regularized(x, a, epsilon, max_iterations)
  else
    n = 0
    an = 1 / a
    sum = an
    while an.abs > epsilon && n < max_iterations
      n += 1
      an *= x / (a + n)
      sum += an
    end
    if n >= max_iterations
      raise Errno::ERANGE
    else
      exp(-x + a * log(x) - log_gamma(a)) * sum
    end
  end
rescue Errno::ERANGE, Errno::EDOM
  0 / 0.0
end

.gammaQ_regularized(x, a, epsilon = 1E-16, max_iterations = 1 << 16) ⇒ Object

Return an approximation of the regularized gammaQ function for x and a with an error of <= epsilon, but only iterate max_iterations-times.



808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
 
# File 'lib/bullshit.rb', line 808

def gammaQ_regularized(x, a, epsilon = 1E-16, max_iterations = 1 << 16)
  x, a = x.to_f, a.to_f
  case
  when a.nan? || x.nan? || a <= 0 || x < 0
    0 / 0.0
  when x == 0
    1.0
  when a > x || a < 1
    1 - gammaP_regularized(x, a, epsilon, max_iterations)
  else
    fraction = ContinuedFraction.for_a do |n, x|
      (2 * n + 1) - a + x
    end.for_b do |n, x|
      n * (a - n)
    end
    exp(-x + a * log(x) - log_gamma(a)) *
      fraction[x, epsilon, max_iterations] ** -1
  end
rescue Errno::ERANGE, Errno::EDOM
  0 / 0.0
end

.log_beta(a, b) ⇒ Object

Returns the natural logarithm of the beta function value for (a, b).



741
742
743
744
745
 
# File 'lib/bullshit.rb', line 741

def log_beta(a, b)
  log_gamma(a) + log_gamma(b) - log_gamma(a + b)
rescue Errno::ERANGE, Errno::EDOM
  0 / 0.0
end

Instance Method Details

#erf(x) ⇒ Object

Returns an approximate value for the error function’s value for x.



835
836
837
838
 
# File 'lib/bullshit.rb', line 835

def erf(x)
  r = sqrt(1 - exp(-x ** 2 * (4 / Math::PI + A * x ** 2) / (1 + A * x ** 2)))
  x < 0 ? -r : r
end

#log_gamma(x) ⇒ Object 



719
720
721
 
# File 'lib/bullshit.rb', line 719

def log_gamma(x)
  lgamma(x).first
end