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.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



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

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



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

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.



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

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



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

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

Raises:

  • (ArgumentError) 


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

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



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

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

Raises:

  • (ArgumentError) 


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

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



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

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(12960) - eps*589.quo(30240))))
    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.



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

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
  return 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



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

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 



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

def q_continued_fraction_complement a,x,with_error=false
  q = q_continued_fraction(a,x,with_error)
  return 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



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

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 



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

def q_large_x_complement a,x,with_error=false
  q = q_large_x(a,x,with_error)
  return 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 



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

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



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

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

Raises:

  • (ArgumentError) 


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

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