1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
|
# Added by John O. Woods, SciRuby project.
module Distribution
module Gamma
module Ruby_
class << self
include Math
# Gamma distribution probability density function
#
# If you're looking at Wikipedia's Gamma distribution page, the arguments for this pdf function correspond
# as follows:
#
# * +x+: same
# * +a+: alpha or k
# + +b+: theta or 1/beta
#
# This is confusing! But we're trying to most closely mirror the GSL function for the gamma distribution
# (see references).
#
# Adapted the function itself from GSL-1.9 in rng/gamma.c: gsl_ran_gamma_pdf
#
# ==References
# * http://www.gnu.org/software/gsl/manual/html_node/The-Gamma-Distribution.html
# * http://en.wikipedia.org/wiki/Gamma_distribution
def pdf(x,a,b)
return 0 if x < 0
if x == 0
return 1.quo(b) if a == 1
return 0
elsif a == 1
Math.exp(-x.quo(b)).quo(b)
else
Math.exp((a-1)*Math.log(x.quo(b)) - x.quo(b) - Math.lgamma(a).first).quo(b)
end
end
# Gamma cumulative distribution function
def cdf(x,a,b)
return 0.0 if x <= 0.0
y = x.quo(b)
return (1-Math::IncompleteGamma.q(a, y)) if y > a
return (Math::IncompleteGamma.p(a, y))
end
#def p_value(pr,a,b)
# cdf(1.0-pr,a,b)
#end
end
end
end
end
|
