Module: Distribution::MathExtension::IncompleteBeta

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

Overview

Calculate regularized incomplete beta function

Constant Summary collapse

MAX_ITER = 
512
CUTOFF = 
2.0 * Float::MIN

Class Method Summary collapse

Class Method Details

.axpy(aa, yy, a, b, x) ⇒ Object

Evaluate aa * beta_inc(a,b,x) + yy

No error mode available.

From GSL-1.9: cdf/beta_inc.c, beta_inc_AXPY



80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
 
# File 'lib/distribution/math_extension/incomplete_beta.rb', line 80

def axpy(aa, yy, a, b, x)
  return aa * 0 + yy if x == 0.0
  return aa * 1 + yy if x == 1.0

  ln_beta   = Math.logbeta(a, b)
  ln_pre    = -ln_beta + a * Math.log(x) + b * Math::Log.log1p(-x)
  prefactor = Math.exp(ln_pre)

  if x < (a + 1).quo(a + b + 2)
    # Apply continued fraction directly
    epsabs  = yy.quo((aa * prefactor).quo(a)).abs * Float::EPSILON
    cf      = continued_fraction(a, b, x, epsabs)
    return  aa * (prefactor * cf).quo(a) + yy
  else
    # Apply continued fraction after hypergeometric transformation
    epsabs = (aa + yy).quo((aa * prefactor).quo(b)) * Float::EPSILON
    cf     = continued_fraction(b, a, 1 - x, epsabs)
    term   = (prefactor * cf).quo(b)
    return aa == -yy ? -aa * term : aa * (1 - term) + yy
  end
end

.continued_fraction(a, b, x, epsabs = nil, with_error = false) ⇒ Object

Continued fraction calculation of incomplete beta beta_cont_frac from GSL-1.9

If epsabs is set, will execute the version of the GSL function in the cdf folder. Otherwise, does the basic one in specfunc.



163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
 
# File 'lib/distribution/math_extension/incomplete_beta.rb', line 163

def continued_fraction(a, b, x, epsabs = nil, with_error = false)
  num_term = 1
  den_term = 1 - (a + b) * x.quo(a + 1)
  k        = 0

  den_term = continued_fraction_cutoff(epsabs) if den_term.abs < CUTOFF
  den_term = 1.quo(den_term)
  cf       = den_term

  1.upto(MAX_ITER) do |k|
    coeff      = k * (b - k) * x.quo(((a - 1) + 2 * k) * (a + 2 * k)) # coefficient for step 1
    delta_frac = nil
    2.times do
      den_term    = 1 + coeff * den_term
      num_term    = 1 + coeff.quo(num_term)

      den_term = continued_fraction_cutoff(epsabs) if den_term.abs < CUTOFF
      num_term = continued_fraction_cutoff(epsabs) if num_term.abs < CUTOFF

      den_term = 1.quo(den_term)

      delta_frac  = den_term * num_term
      cf *= delta_frac

      coeff = -(a + k) * (a + b + k) * x.quo((a + 2 * k) * (a + 2 * k + 1)) # coefficient for step 2
    end

    break if (delta_frac - 1).abs < 2.0 * Float::EPSILON
    break if !epsabs.nil? && (cf * (delta_frac - 1).abs < epsabs)
  end

  if k > MAX_ITER
    fail('Exceeded maximum number of iterations') if epsabs.nil?
    return with_error ? [0.0 / 0, 0] : 0.0 / 0 # NaN if epsabs is set
  end

  with_error ? [cf, k * 4 * Float::EPSILON * cf.abs] : cf
end

.continued_fraction_cutoff(epsabs) ⇒ Object 



153
154
155
156
 
# File 'lib/distribution/math_extension/incomplete_beta.rb', line 153

def continued_fraction_cutoff(epsabs)
  return CUTOFF if epsabs.nil?
  0.0 / 0 # NaN
end

.evaluate(a, b, x, with_error = false) ⇒ Object

Evaluate the incomplete beta function gsl_sf_beta_inc_e



104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
 
# File 'lib/distribution/math_extension/incomplete_beta.rb', line 104

def evaluate(a, b, x, with_error = false)
  fail(ArgumentError, "Domain error: a(#{a}), b(#{b}) must be positive; x(#{x}) must be between 0 and 1, inclusive") if a <= 0 || b <= 0 || x < 0 || x > 1
  if x == 0
    return with_error ? [0.0, 0.0] : 0.0
  elsif x == 1
    return with_error ? [1.0, 0.0] : 1.0
  else

    ln_beta = Beta.log_beta(a, b, with_error)
    ln_1mx  = Log.log_1plusx(-x, with_error)
    ln_x    = Math.log(x)

    ln_beta, ln_beta_err, ln_1mx, ln_1mx_err, ln_x_err = begin
      # STDERR.puts("Warning: Error is unknown for Math::log, guessing.")
      [ln_beta, ln_1mx, Float::EPSILON].flatten
    end

    ln_pre      = -ln_beta + a * ln_x + b * ln_1mx
    ln_pre_err  = ln_beta_err + (a * ln_x_err).abs + (b * ln_1mx_err).abs if with_error

    prefactor, prefactor_err   = begin
      if with_error
        exp_err(ln_pre, ln_pre_err)
      else
        [Math.exp(ln_pre), nil]
      end
    end

    if x < (a + 1).quo(a + b + 2)
      # Apply continued fraction directly

      cf      = continued_fraction(a, b, x, nil, with_error)
      cf, cf_err = cf if with_error
      result  = (prefactor * cf).quo(a)
      return with_error ? [result, ((prefactor_err * cf).abs + (prefactor * cf_err).abs).quo(a)] : result
    else
      # Apply continued fraction after hypergeometric transformation

      cf      = continued_fraction(b, a, 1 - x, nil)
      cf, cf_err = cf if with_error
      term    = (prefactor * cf).quo(b)
      result  = 1 - term

      return with_error ? [result, (prefactor_err * cf).quo(b) + (prefactor * cf_err).quo(b) + 2.0 * Float::EPSILON * (1 + term.abs)] : result
    end

  end
end