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



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

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.



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
201
202
203
204
205
206
 
# File 'lib/distribution/math_extension/incomplete_beta.rb', line 168

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



158
159
160
161
 
# File 'lib/distribution/math_extension/incomplete_beta.rb', line 158

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

Raises:

  • (ArgumentError) 


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
152
153
154
155
 
# File 'lib/distribution/math_extension/incomplete_beta.rb', line 108

def evaluate(a,b,x,with_error=false)
  raise(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