From 97ae69b60cf9b28d524f07b1effeb29c872cb734 Mon Sep 17 00:00:00 2001
From: Simon Byrne <simonbyrne@gmail.com>
Date: Wed, 19 Nov 2014 15:49:31 +0000
Subject: [PATCH] reconfigure samplers, new gamma sampler

---
 src/Distributions.jl  |   2 +-
 src/samplers/gamma.jl | 137 +++++++++++++++++++++++++++++-------------
 test/samplers.jl      |  23 +++++--
 3 files changed, 112 insertions(+), 50 deletions(-)

diff --git a/src/Distributions.jl b/src/Distributions.jl
index 662372f30..1473b07b1 100644
--- a/src/Distributions.jl
+++ b/src/Distributions.jl
@@ -243,11 +243,11 @@ include("functionals.jl")
 include("genericfit.jl")
 
 # specific samplers and distributions
-include("samplers.jl")
 include("univariates.jl")
 include("empirical.jl")
 include("multivariates.jl")
 include("matrixvariates.jl")
+include("samplers.jl")
 
 # others
 include("truncate.jl")
diff --git a/src/samplers/gamma.jl b/src/samplers/gamma.jl
index becd476ee..875260d4f 100644
--- a/src/samplers/gamma.jl
+++ b/src/samplers/gamma.jl
@@ -1,12 +1,10 @@
 
 immutable GammaRmathSampler <: Sampleable{Univariate,Continuous}
-    α::Float64
-    scale::Float64
+    d::Gamma
 end
 
-rand(s::GammaRmathSampler) = 
-    ccall((:rgamma, "libRmath-julia"), Float64, (Float64, Float64), s.α, s.scale)
-
+rand(s::GammaRmathSampler) =
+    ccall((:rgamma, "libRmath-julia"), Float64, (Float64, Float64), s.d.shape, s.d.scale)
 
 
 # "Generating gamma variates by a modified rejection technique"
@@ -14,12 +12,13 @@ rand(s::GammaRmathSampler) =
 # Communications of the ACM, Vol 25(1), 1982, pp 47-54
 # doi:10.1145/358315.358390
 
-# suitable for scale >= 1.0
+# suitable for shape >= 1.0
 
 immutable GammaGDSampler <: Sampleable{Univariate,Continuous}
     a::Float64
     s2::Float64
     s::Float64
+    i2s::Float64
     d::Float64
     q0::Float64
     b::Float64
@@ -28,10 +27,13 @@ immutable GammaGDSampler <: Sampleable{Univariate,Continuous}
     scale::Float64
 end
 
-function GammaGDSampler(a::Float64,scale::Float64)
+function GammaGDSampler(g::Gamma)
+    a = g.shape
+
     # Step 1
     s2 = a-0.5
     s = sqrt(s2)
+    i2s = 0.5/s
     d = 5.656854249492381 - 12.0s # 4*sqrt(2) - 12s
 
     # Step 4
@@ -61,7 +63,7 @@ function GammaGDSampler(a::Float64,scale::Float64)
         c = 0.1515/s
     end
 
-    GammaGDSampler(a,s2,s,d,q0,b,σ,c,scale)
+    GammaGDSampler(a,s2,s,i2s,d,q0,b,σ,c,g.scale)
 end
 
 function rand(s::GammaGDSampler)
@@ -77,7 +79,7 @@ function rand(s::GammaGDSampler)
     # Step 5
     if x > 0.0
         # Step 6
-        v = t/(2.0*s.s)
+        v = t*s.i2s
         if abs(v) > 0.25
             q = s.q0 - s.s*t + 0.25*t*t + 2.0*s.s2*log1p(v)
         else
@@ -107,7 +109,7 @@ function rand(s::GammaGDSampler)
     t < -0.718_744_837_717_19 && @goto step8
 
     # Step 10
-    v = t/(2.0*s.s)
+    v = t*s.i2s
     if abs(v) > 0.25
         q = s.q0 - s.s*t + 0.25*t*t + 2.0*s.s2*log1p(v)
     else
@@ -131,61 +133,110 @@ function rand(s::GammaGDSampler)
     return x*x*s.scale
 end
 
-#  A simple method for generating gamma variables - Marsaglia and Tsang (2000)
-#  http://www.cparity.com/projects/AcmClassification/samples/358414.pdf
-#  Page 369
-#  basic simulation loop for pre-computed d and c
-#
+# "Computer methods for sampling from gamma, beta, poisson and bionomial distributions"
+# J.H. Ahrens and U. Dieter
+# Computing, 1974, Volume 12(3), pp 223-246
+# doi:10.1007/BF02293108
+
+# valid for 0 < shape <= 1
+immutable GammaGSSampler <: Sampleable{Univariate,Continuous}
+    a::Float64
+    ia::Float64
+    b::Float64
+    scale::Float64
+end
+
+function GammaGSSampler(d::Gamma)
+    a = d.shape
+    ia = 1/d.shape
+    b = 1.0+0.36787944117144233*d.shape
+    GammaGSSampler(a,ia,b,d.scale)
+end
+
+function rand(s::GammaGSSampler)
+    while true
+        # step 1
+        p = s.b*rand()
+        e = Base.Random.randmtzig_exprnd()
+        if p <= 1.0
+            # step 2
+            x = exp(log(p)*s.ia)
+            e < x || return s.scale*x
+        else
+            # step 3
+            x = -log(s.ia*(s.b-p))
+            e < log(x)*(1.0-s.a) || return s.scale*x
+        end
+    end
+end
+
+
+# "A simple method for generating gamma variables"
+# G. Marsaglia and W.W. Tsang
+# ACM Transactions on Mathematical Software (TOMS), 2000, Volume 26(3), pp. 363-372
+# doi:10.1145/358407.358414
+# http://www.cparity.com/projects/AcmClassification/samples/358414.pdf
 
 immutable GammaMTSampler <: Sampleable{Univariate,Continuous}
-    iα::Float64
     d::Float64
     c::Float64
     κ::Float64
 end
 
-function GammaMTSampler(α::Float64, scale::Float64)
-    if α >= 1.0
-        iα = 1.0
-        d = α - 1/3
-    else
-        iα = 1.0 / α
-        d = α + 2/3
-    end
+function GammaMTSampler(g::Gamma)
+    d = g.shape - 1/3
     c = 1.0 / sqrt(9.0 * d)
-    κ = d * scale
-    GammaMTSampler(iα, d, c, κ)
+    κ = d * g.scale
+    GammaMTSampler(d, c, κ)
 end
 
-GammaMTSampler(α::Float64) = GammaMTSampler(α, 1.0)
-
 function rand(s::GammaMTSampler)
-    d = s.d
-    c = s.c
-    iα = s.iα
-
-    v = 0.0
     while true
         x = randn()
-        v = 1.0 + c * x
+        v = 1.0 + s.c * x
         while v <= 0.0
             x = randn()
-            v = 1.0 + c * x
+            v = 1.0 + s.c * x
         end
         v *= (v * v)
         u = rand()
         x2 = x * x
         if u < 1.0 - 0.331 * abs2(x2)
-            break
+            return v*s.κ
         end
-        if log(u) < 0.5 * x2 + d * (1.0 - v + log(v))
-            break
+        if log(u) < 0.5 * x2 + s.d * (1.0 - v + log(v)) # logmxp1
+            return v*s.κ
         end
     end
-    v *= s.κ
-    if iα > 1.0
-        v *= (rand() ^ iα)
-    end
-    return v
 end
 
+# Inverse Power sampler
+# uses the x*u^(1/a) trick from Marsaglia and Tsang (2000) for when shape < 1
+immutable GammaIPSampler{S<:Sampleable{Univariate,Continuous}} <: Sampleable{Univariate,Continuous}
+    s::S #sampler for Gamma(1+shape,scale)
+    nia::Float64 #-1/scale
+end
+
+function GammaIPSampler{S<:Sampleable}(d::Gamma,::Type{S})
+    GammaIPSampler(Gamma(1.0+d.shape,d.scale), -1.0/d.shape)
+end
+GammaIPSampler(d::Gamma) = GammaIPSampler(d,GammaGDSampler)
+
+function rand(s::GammaIPSampler)
+    x = rand(s.s)
+    e = Base.Random.randmtzig_exprnd()
+    x*exp(s.nia*e)
+end
+
+# function sampler(d::Gamma)
+#     if d.shape < 1.0
+#         # TODO: d.shape = 0.5 : use scaled chisq
+#         GammaIPSampler(d,GammaGDSampler)
+#     elseif d.shape == 1.0
+#         Exponential(d.scale)
+#     else
+#         GammaGDSampler(d)
+#     end
+# end
+
+# rand(d::Gamma) = rand(sampler(d))
diff --git a/test/samplers.jl b/test/samplers.jl
index 649c3916a..b5462fa8a 100644
--- a/test/samplers.jl
+++ b/test/samplers.jl
@@ -13,7 +13,10 @@ import Distributions:
     PoissonADSampler, 
     PoissonCountSampler,
     ExponentialSampler,
-    GammaMTSampler
+    GammaGDSampler,
+    GammaGSSampler,
+    GammaMTSampler,
+    GammaIPSampler
 
 n_tsamples = 10^6
 
@@ -74,11 +77,19 @@ end
 
 
 ## Gamma samplers
-
-for S in [GammaMTSampler]
-    for pa in [(1.0, 1.0), (2.0, 1.0), (3.0, 1.0), (0.5, 1.0), 
-               (1.0, 2.0), (3.0, 2.0), (0.5, 2.0)]
-        test_samples(S(pa...), Gamma(pa...), n_tsamples)
+# shape >= 1
+for S in [GammaGDSampler, GammaMTSampler]
+    println("    testing $S")
+    for d in [Gamma(1.0, 1.0), Gamma(2.0, 1.0), Gamma(3.0, 1.0),
+               Gamma(1.0, 2.0), Gamma(3.0, 2.0), Gamma(100.0, 2.0)]
+        test_samples(S(d), d, n_tsamples)
     end
 end
 
+# shape < 1
+for S in [GammaGSSampler, GammaIPSampler]
+    println("    testing $S")
+    for d in [Gamma(0.1,1.0),Gamma(0.9,1.0)]
+        test_samples(S(d), d, n_tsamples)
+    end
+end