# 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

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

• Continued fraction calculation of incomplete beta beta_cont_frac from GSL-1.9.

• Evaluate the incomplete beta function gsl_sf_beta_inc_e.

## 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```