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.

Raises:

  • (ArgumentError) 


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
67
 
# File 'lib/distribution/math_extension/incomplete_beta.rb', line 10

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

  raise(ArgumentError, "x and y must be nonzero") if x == 0.0 || y == 0.0
  raise(ArgumentError, "not defined for negative integers") if [x,y].any? { |v| (v.is_a?(Fixnum) && 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

  raise("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