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
|
# Added by John O. Woods, SciRuby project.
module Distribution
module Beta
module Ruby_
class << self
include Math
# Beta distribution probability density function
#
# Adapted from GSL-1.9 (apparently by Knuth originally), found in randist/beta.c
#
# Form: p(x) dx = (Gamma(a + b)/(Gamma(a) Gamma(b))) x^(a-1) (1-x)^(b-1) dx
#
# == References
# * http://www.gnu.org/s/gsl/manual/html_node/The-Gamma-Distribution.html
def pdf(x,a,b)
return 0 if x < 0 || x > 1
gab = Math.lgamma(a+b).first
ga = Math.lgamma(a).first
gb = Math.lgamma(b).first
if x == 0.0 || x == 1.0
Math.exp(gab - ga - gb) * x**(a-1) * (1-x)**(b-1)
else
Math.exp(gab - ga - gb + Math.log(x)*(a-1) + Math::Log.log1p(-x)*(b-1))
end
end
# Gamma cumulative distribution function
# Translated from GSL-1.9: cdf/beta.c gsl_cdf_beta_P
def cdf(x,a,b)
return 0.0 if x <= 0.0
return 1.0 if x >= 1.0
Math::IncompleteBeta.axpy(1.0, 0.0, a,b,x)
end
end
end
end
end
|
