Module: Distribution::MathExtension::IncompleteGamma

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

Constant Summary collapse

NMAX = 
5000
SMALL = 
Float::EPSILON**3
PG21 =

PolyGamma[2,1]

-2.404113806319188570799476 # PolyGamma[2,1]

Class Method Summary collapse

Class Method Details

.a_greater_than_0(a, x, with_error = false) ⇒ Object

gamma_inc_a_gt_0



355
356
357
358
359
360
361
362
363
364
 
# File 'lib/distribution/math_extension/incomplete_gamma.rb', line 355

def a_greater_than_0(a, x, with_error = false)
  q       = q(a, x, with_error)
  q, q_err = q if with_error
  g       = Math.gamma(a)
  STDERR.puts("Warning: Don't know error for Math.gamma. Error will be incorrect") if with_error
  g_err   = Float::EPSILON
  result  = g * q
  error   = (g * q_err).abs + (g_err * q).abs if with_error
  with_error ? [result, error] : result
end

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

gamma_inc_CF



367
368
369
370
371
372
373
374
375
376
377
378
379
380
 
# File 'lib/distribution/math_extension/incomplete_gamma.rb', line 367

def continued_fraction(a, x, with_error = false)
  f = f_continued_fraction(a, x, with_error)
  f, f_error = f if with_error
  pre = Math.exp((a - 1.0) * Math.log(x) - x)
  STDERR.puts("Warning: Don't know error for Math.exp. Error will be incorrect") if with_error
  pre_error = Float::EPSILON
  result    = f * pre
  if with_error
    error     = (f_error * pre).abs + (f * pre_error) + (2.0 + a.abs) * Float::EPSILON * result.abs
    [result, error]
  else
    result
  end
end

.d(a, x, with_error = false) ⇒ Object

The dominant part, D(a,x) := x^a e^(-x) / Gamma(a+1) gamma_inc_D in GSL-1.9.



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
 
# File 'lib/distribution/math_extension/incomplete_gamma.rb', line 38

def d(a, x, with_error = false)
  error = nil
  if a < 10.0
    ln_a = Math.lgamma(a + 1.0).first
    lnr  = a * Math.log(x) - x - ln_a
    result = Math.exp(lnr)
    error = 2.0 * Float::EPSILON * (lnr.abs + 1.0) + result.abs if with_error
    with_error ? [result, error] : result
  else
    ln_term = ln_term_error = nil
    if x < 0.5 * a
      u       = x / a.to_f
      ln_u    = Math.log(u)
      ln_term = ln_u - u + 1.0
      ln_term_error = (ln_u.abs + u.abs + 1.0) * Float::EPSILON if with_error
    else
      mu      = (x - a) / a.to_f
      ln_term = Log.log_1plusx_minusx(mu, with_error)
      ln_term, ln_term_error = ln_term if with_error
    end
    gstar = Gammastar.evaluate(a, with_error)
    gstar, gstar_error = gstar if with_error
    term1 = Math.exp(a * ln_term) / Math.sqrt(2.0 * Math::PI * a)
    result = term1 / gstar
    error  = 2.0 * Float::EPSILON * ((a * ln_term).abs + 1.0) * result.abs + gstar_error / gstar.abs * result.abs if with_error
    with_error ? [result, error] : result
  end
end

.f_continued_fraction(a, x, with_error = false) ⇒ Object

gamma_inc_F_CF



231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
 
# File 'lib/distribution/math_extension/incomplete_gamma.rb', line 231

def f_continued_fraction(a, x, with_error = false)
  hn = 1.0 # convergent
  cn = 1.0 / SMALL
  dn = 1.0
  n  = 2
  2.upto(NMAX - 1).each do |n|
    an = n.odd? ? 0.5 * (n - 1) / x : (0.5 * n - a) / x
    dn = 1.0 + an * dn
    dn = SMALL if dn.abs < SMALL
    cn = 1.0 + an / cn
    cn = SMALL if cn.abs < SMALL
    dn = 1.0 / dn
    delta = cn * dn
    hn *= delta
    break if (delta - 1.0).abs < Float::EPSILON
  end

  if n == NMAX
    STDERR.puts('Error: n reached NMAX in f continued fraction')
  else
    with_error ? [hn, 2.0 * Float::EPSILON * hn.abs + Float::EPSILON * (2.0 + 0.5 * n) * hn.abs] : hn
  end
end

.p(a, x, with_error = false) ⇒ Object

The incomplete gamma function. gsl_sf_gamma_inc_P_e



109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
 
# File 'lib/distribution/math_extension/incomplete_gamma.rb', line 109

def p(a, x, with_error = false)
  fail(ArgumentError, 'Range Error: a must be positive, x must be non-negative') if a <= 0.0 || x < 0.0
  if x == 0.0
    return with_error ? [0.0, 0.0] : 0.0
  elsif x < 20.0 || x < 0.5 * a
    return p_series(a, x, with_error)
  elsif a > 1e6 && (x - a) * (x - a) < a
    return q_asymptotic_uniform_complement a, x, with_error
  elsif a <= x
    if a > 0.2 * x
      return q_continued_fraction_complement(a, x, with_error)
    else
      return q_large_x_complement(a, x, with_error)
    end
  elsif (x - a) * (x - a) < a
    return q_asymptotic_uniform_complement a, x, with_error
  else
    return p_series(a, x, with_error)
  end
end

.p_series(a, x, with_error = false) ⇒ Object

gamma_inc_P_series



68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
 
# File 'lib/distribution/math_extension/incomplete_gamma.rb', line 68

def p_series(a, x, with_error = false)
  d = d(a, x, with_error)
  d, d_err = d if with_error
  sum      = 1.0
  term     = 1.0
  n        = 1
  1.upto(NMAX - 1) do |n|
    term *= x / (a + n).to_f
    sum += term
    break if (term / sum).abs < Float::EPSILON
  end

  result   = d * sum

  if n == NMAX
    STDERR.puts('Error: n reached NMAX in p series')
  else
    return with_error ? [result, d_err * sum.abs + (1.0 + n) * Float::EPSILON * result.abs] : result
  end
end

.q(a, x, with_error = false) ⇒ Object

gamma_inc_Q_e



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
156
157
158
159
 
# File 'lib/distribution/math_extension/incomplete_gamma.rb', line 131

def q(a, x, with_error = false)
  fail(ArgumentError, 'Range Error: a and x must be non-negative') if a < 0.0 || x < 0.0
  if x == 0.0
    return with_error ? [1.0, 0.0] : 1.0
  elsif a == 0.0
    return with_error ? [0.0, 0.0] : 0.0
  elsif x <= 0.5 * a
    # If series is quick, do that.
    p = p_series(a, x, with_error)
    p, p_err = p if with_error
    result  = 1.0 - p
    return with_error ? [result, p_err + 2.0 * Float::EPSILON * result.abs] : result
  elsif a >= 1.0e+06 && (x - a) * (x - a) < a # difficult asymptotic regime, only way to do this region
    return q_asymptotic_uniform(a, x, with_error)
  elsif a < 0.2 && x < 5.0
    return q_series(a, x, with_error)
  elsif a <= x
    return x <= 1.0e+06 ? q_continued_fraction(a, x, with_error) : q_large_x(a, x, with_error)
  else
    if x > a - Math.sqrt(a)
      return q_continued_fraction(a, x, with_error)
    else
      p = p_series(a, x, with_error)
      p, p_err = p if with_error
      result = 1.0 - p
      return with_error ? [result, p_err + 2.0 * Float::EPSILON * result.abs] : result
    end
  end
end

.q_asymptotic_uniform(a, x, with_error = false) ⇒ Object

Uniform asymptotic for x near a, a and x large gamma_inc_Q_asymp_unif



201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
 
# File 'lib/distribution/math_extension/incomplete_gamma.rb', line 201

def q_asymptotic_uniform(a, x, with_error = false)
  rta = Math.sqrt(a)
  eps = (x - a).quo(a)

  ln_term = Log.log_1plusx_minusx(eps, with_error)
  ln_term, ln_term_err = ln_term if with_error

  eta     = (eps >= 0 ? 1 : -1) * Math.sqrt(-2 * ln_term)

  erfc    = Math.erfc_e(eta * rta / SQRT2, with_error)
  erfc, erfc_err = erfc if with_error

  c0 = c1 = nil
  if eps.abs < ROOT5_FLOAT_EPSILON
    c0 = -1.quo(3) + eps * (1.quo(12) - eps * (23.quo(540) - eps * (353.quo(12_960) - eps * 589.quo(30_240))))
    c1 = -1.quo(540) - eps.quo(288)
  else
    rt_term = Math.sqrt(-2 * ln_term.quo(eps * eps))
    lam     = x.quo(a)
    c0      = (1 - 1 / rt_term) / eps
    c1      = -(eta**3 * (lam * lam + 10 * lam + 1) - 12 * eps**3).quo(12 * eta**3 * eps**3)
  end

  r = Math.exp(-0.5 * a * eta * eta) / (SQRT2 * SQRTPI * rta) * (c0 + c1.quo(a))

  result = 0.5 * erfc + r
  with_error ? [result, Float::EPSILON + (r * 0.5 * a * eta * eta).abs + 0.5 * erfc_err + 2.0 * Float::EPSILON + result.abs] : result
end

.q_asymptotic_uniform_complement(a, x, with_error = false) ⇒ Object

This function does not exist in GSL, but is nonetheless GSL code. It's for calculating two specific ranges of p.



90
91
92
93
94
95
 
# File 'lib/distribution/math_extension/incomplete_gamma.rb', line 90

def q_asymptotic_uniform_complement(a, x, with_error = false)
  q = q_asymptotic_uniform(a, x, with_error)
  q, q_err = q if with_error
  result = 1.0 - q
  with_error ? [result, q_err + 2.0 * Float::EPSILON * result.abs] : result
end

.q_continued_fraction(a, x, with_error = false) ⇒ Object

gamma_inc_Q_CF



162
163
164
165
166
167
168
169
170
171
 
# File 'lib/distribution/math_extension/incomplete_gamma.rb', line 162

def q_continued_fraction(a, x, with_error = false)
  d = d(a, x, with_error)
  f = f_continued_fraction(a, x, with_error)

  if with_error
    [d.first * (a / x).to_f * f.first, d.last * ((a / x).to_f * f.first).abs + (d.first * a / x * f.last).abs]
  else
    d * (a / x).to_f * f
  end
end

.q_continued_fraction_complement(a, x, with_error = false) ⇒ Object 



97
98
99
100
 
# File 'lib/distribution/math_extension/incomplete_gamma.rb', line 97

def q_continued_fraction_complement(a, x, with_error = false)
  q = q_continued_fraction(a, x, with_error)
  with_error ? [1.0 - q.first, q.last + 2.0 * Float::EPSILON * (1.0 - q.first).abs] : 1.0 - q
end

.q_large_x(a, x, with_error = false) ⇒ Object

gamma_inc_Q_large_x in GSL-1.9



174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
 
# File 'lib/distribution/math_extension/incomplete_gamma.rb', line 174

def q_large_x(a, x, with_error = false)
  d = d(a, x, with_error)
  d, d_err = d if with_error
  sum  = 1.0
  term = 1.0
  last = 1.0
  n    = 1
  1.upto(NMAX - 1).each do |n|
    term *= (a - n) / x
    break if (term / last).abs > 1.0
    break if (term / sum).abs < Float::EPSILON
    sum += term
    last  = term
  end

  result = d * (a / x) * sum
  error  = d_err * (a / x).abs * sum if with_error

  if n == NMAX
    STDERR.puts('Error: n reached NMAX in q_large_x')
  else
    return with_error ? [result, error] : result
  end
end

.q_large_x_complement(a, x, with_error = false) ⇒ Object 



102
103
104
105
 
# File 'lib/distribution/math_extension/incomplete_gamma.rb', line 102

def q_large_x_complement(a, x, with_error = false)
  q = q_large_x(a, x, with_error)
  with_error ? [1.0 - q.first, q.last + 2.0 * Float::EPSILON * (1.0 - q.first).abs] : 1.0 - q
end

.q_series(a, x, with_error = false) ⇒ Object 



255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
 
# File 'lib/distribution/math_extension/incomplete_gamma.rb', line 255

def q_series(a, x, with_error = false)
  term1 = nil
  sum   = nil
  term2 = nil
  begin
    lnx  = Math.log(x)
    el   = EULER + lnx
    c1   = -el
    c2   = Math::PI * Math::PI / 12.0 - 0.5 * el * el
    c3   = el * (Math::PI * Math::PI / 12.0 - el * el / 6.0) + PG21 / 6.0
    c4   = -0.04166666666666666667 *
           (-1.758243446661483480 + lnx) *
           (-0.764428657272716373 + lnx) *
           (0.723980571623507657 + lnx) *
           (4.107554191916823640 + lnx)
    c5 = -0.0083333333333333333 *
         (-2.06563396085715900 + lnx) *
         (-1.28459889470864700 + lnx) *
         (-0.27583535756454143 + lnx) *
         (1.33677371336239618 + lnx) *
         (5.17537282427561550 + lnx)
    c6 = -0.0013888888888888889 *
         (-2.30814336454783200 + lnx) *
         (-1.65846557706987300 + lnx) *
         (-0.88768082560020400 + lnx) *
         (0.17043847751371778 + lnx) *
         (1.92135970115863890 + lnx) *
         (6.22578557795474900 + lnx)
    c7 = -0.00019841269841269841
    (-2.5078657901291800 + lnx) *
      (-1.9478900888958200 + lnx) *
      (-1.3194837322612730 + lnx) *
      (-0.5281322700249279 + lnx) *
      (0.5913834939078759 + lnx) *
      (2.4876819633378140 + lnx) *
      (7.2648160783762400 + lnx)
    c8 = -0.00002480158730158730 *
         (-2.677341544966400 + lnx) *
         (-2.182810448271700 + lnx) *
         (-1.649350342277400 + lnx) *
         (-1.014099048290790 + lnx) *
         (-0.191366955370652 + lnx) *
         (0.995403817918724 + lnx) *
         (3.041323283529310 + lnx) *
         (8.295966556941250 + lnx) *
         c9 = -2.75573192239859e-6 *
              (-2.8243487670469080 + lnx) *
              (-2.3798494322701120 + lnx) *
              (-1.9143674728689960 + lnx) *
              (-1.3814529102920370 + lnx) *
              (-0.7294312810261694 + lnx) *
              (0.1299079285269565 + lnx) *
              (1.3873333251885240 + lnx) *
              (3.5857258865210760 + lnx) *
              (9.3214237073814600 + lnx) *
              c10 = -2.75573192239859e-7 *
                    (-2.9540329644556910 + lnx) *
                    (-2.5491366926991850 + lnx) *
                    (-2.1348279229279880 + lnx) *
                    (-1.6741881076349450 + lnx) *
                    (-1.1325949616098420 + lnx) *
                    (-0.4590034650618494 + lnx) *
                    (0.4399352987435699 + lnx) *
                    (1.7702236517651670 + lnx) *
                    (4.1231539047474080 + lnx) *
                    (10.342627908148680 + lnx)
    term1 = a * (c1 + a * (c2 + a * (c3 + a * (c4 + a * (c5 + a * (c6 + a * (c7 + a * (c8 + a * (c9 + a * c10)))))))))
  end

  n   = 1
  begin
    t   = 1.0
    sum = 1.0
    1.upto(NMAX - 1).each do |n|
      t *= -x / (n + 1.0)
      sum += (a + 1.0) / (a + n + 1.0) * t
      break if (t / sum).abs < Float::EPSILON
    end
  end

  if n == NMAX
    STDERR.puts('Error: n reached NMAX in q_series')
  else
    term2 = (1.0 - term1) * a / (a + 1.0) * x * sum
    result = term1 + term2
    with_error ? [result, Float::EPSILON * term1.abs + 2.0 * term2.abs + 2.0 * Float::EPSILON * result.abs] : result
  end
end

.series(a, x, with_error = false) ⇒ Object

gamma_inc_series



345
346
347
348
349
350
351
352
 
# File 'lib/distribution/math_extension/incomplete_gamma.rb', line 345

def series(a, x, with_error = false)
  q = q_series(a, x, with_error)
  g = Math.gamma(a)
  STDERR.puts("Warning: Don't know error for Math.gamma. Error will be incorrect") if with_error
  # When we get the error from Gamma, switch the comment on the next to lines
  # with_error ? [q.first*g.first, (q.first*g.last).abs + (q.last*g.first).abs + 2.0*Float::EPSILON*(q.first*g.first).abs] : q*g
  with_error ? [q.first * g, (q.first * Float::EPSILON).abs + (q.last * g.first).abs + 2.0 * Float::EPSILON(q.first * g).abs] : q * g
end

.unnormalized(a, x, with_error = false) ⇒ Object

Unnormalized incomplete gamma function. gsl_sf_gamma_inc_e



384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
 
# File 'lib/distribution/math_extension/incomplete_gamma.rb', line 384

def unnormalized(a, x, with_error = false)
  fail(ArgumentError, 'x cannot be negative') if x < 0.0

  if x == 0.0
    result  = Math.gamma(a.to_f)
    STDERR.puts("Warning: Don't know error for Math.gamma. Error will be incorrect") if with_error
    return with_error ? [result, Float::EPSILON] : result
  elsif a == 0.0
    return ExponentialIntegral.first_order(x.to_f, with_error)
  elsif a > 0.0
    return a_greater_than_0(a.to_f, x.to_f, with_error)
  elsif x > 0.25
    # continued fraction seems to fail for x too small
    return continued_fraction(a.to_f, x.to_f, with_error)
  elsif a.abs < 0.5
    return series(a.to_f, x.to_f, with_error)
  else
    fa = a.floor.to_f
    da = a - fa
    g_da = da > 0.0 ? a_greater_than_0(da, x.to_f, with_error) : ExponentialIntegral.first_order(x.to_f, with_error)
    g_da, g_da_err = g_da if with_error
    alpha = da
    gax = g_da

    # Gamma(alpha-1,x) = 1/(alpha-1) (Gamma(a,x) - x^(alpha-1) e^-x)
    begin
      shift  = Math.exp(-x + (alpha - 1.0) * Math.log(x))
      gax    = (gax - shift) / (alpha - 1.0)
      alpha -= 1.0
    end while alpha > a

    result = gax
    return with_error ? [result, 2.0 * (1.0 + a.abs) * Float::EPSILON * gax.abs] : result
  end
end