Module: Distribution::MathExtension::Beta

Defined in:
lib/distribution/math_extension/incomplete_beta.rb

Class Method Summary collapse

Class Method Details

.log_beta(x, y, with_error = false) ⇒ Object

Based on gsl_sf_lnbeta_e and gsl_sf_lnbeta_sgn_e Returns result and sign in an array. If with_error is specified, also returns the error.



10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
 
# File 'lib/distribution/math_extension/incomplete_beta.rb', line 10

def log_beta(x, y, with_error = false)
  sign = nil

  fail(ArgumentError, 'x and y must be nonzero') if x == 0.0 || y == 0.0
  fail(ArgumentError, 'not defined for negative integers') if [x, y].any? { |v| v < 0 }

  # See if we can handle the positive case with min/max < 0.2
  if x > 0 && y > 0
    min, max = [x, y].minmax
    ratio    = min.quo(max)

    if ratio < 0.2
      gsx   = Gammastar.evaluate(x, with_error)
      gsy   = Gammastar.evaluate(y, with_error)
      gsxy  = Gammastar.evaluate(x + y, with_error)
      lnopr = Log.log_1plusx(ratio, with_error)

      gsx, gsx_err, gsy, gsy_err, gsxy, gsxy_err, lnopr, lnopr_err = [gsx, gsy, gsxy, lnopr].flatten if with_error

      lnpre = Math.log((gsx * gsy).quo(gsxy) * Math::SQRT2 * Math::SQRTPI)
      lnpre_err = gsx_err.quo(gsx) + gsy_err(gsy) + gsxy_err.quo(gsxy) if with_error

      t1    = min * Math.log(ratio)
      t2    = 0.5 * Math.log(min)
      t3    = (x + y - 0.5) * lnopr

      lnpow       = t1 - t2 - t3
      lnpow_err   = Float::EPSILON * (t1.abs + t2.abs + t3.abs) + (x + y - 0.5).abs * lnopr_err if with_error

      result      = lnpre + lnpow
      error       = lnpre_err + lnpow_err + 2.0 * Float::EPSILON * result.abs if with_error

      return with_error ? [result, 1.0, error] : [result, 1.0]
    end
  end

  # General case: fallback
  lgx, sgx   = Math.lgamma(x)
  lgy, sgy   = Math.lgamma(y)
  lgxy, sgxy = Math.lgamma(x + y)
  sgn        = sgx * sgy * sgxy

  fail('Domain error: sign is -') if sgn == -1

  result = lgx + lgy - lgxy
  if with_error
    lgx_err, lgy_err, lgxy_err = begin
      STDERR.puts('Warning: Error is unknown for Math::lgamma, guessing.')
      [Math::EPSILON, Math::EPSILON, Math::EPSILON]
    end

    error  = lgx_err + lgy_err + lgxy_err + Float::EPSILON * (lgx.abs + lgy.abs + lgxy.abs) + 2.0 * (Float::EPSILON) * result.abs
    return [result, sgn, error]
  else
    return [result, sgn]
  end
end