Merge 685b7e3 into a314e5c

.travis.yml

@@ -1,8 +1,10 @@
 language: julia
 julia:
-  - 0.6.0
-  - 0.6.1
-  - 0.6.2
+  - 0.7.0
+  - 1.0.0
+  - 1.0.1
+  - 1.0.2
+  - 1.0.3
 before_install:
   - pip install --user codecov
 after_success:

README.md

@@ -1,6 +1,6 @@
 [![Build Status](https://travis-ci.org/mauriciogtec/AdaptiveRejectionSampling.jl.svg?branch=master)](https://travis-ci.org/mauriciogtec/AdaptiveRejectionSampling.jl)
-[![AdaptiveRejectionSampling](http://pkg.julialang.org/badges/AdaptiveRejectionSampling_0.6.svg)](http://pkg.julialang.org/detail/AdaptiveRejectionSampling)
 [![AdaptiveRejectionSampling](http://pkg.julialang.org/badges/AdaptiveRejectionSampling_0.7.svg)](http://pkg.julialang.org/detail/AdaptiveRejectionSampling)
+[![AdaptiveRejectionSampling](http://pkg.julialang.org/badges/AdaptiveRejectionSampling_1.0.svg)](http://pkg.julialang.org/detail/AdaptiveRejectionSampling)
 
 [![Coverage Status](https://coveralls.io/repos/github/mauriciogtec/AdaptiveRejectionSampling.jl/badge.svg?branch=master)](https://coveralls.io/github/mauriciogtec/AdaptiveRejectionSampling.jl?branch=master)
 [![codecov](https://codecov.io/gh/mauriciogtec/AdaptiveRejectionSampling.jl/branch/master/graph/badge.svg)](https://codecov.io/gh/mauriciogtec/AdaptiveRejectionSampling.jl)
@@ -36,9 +36,10 @@ Let's verify the result
 
 
 ```julia
-x = -linspace(-10.0, 10.0, 100)
-envelop = eval_envelop.(sampler.envelop, x)
-target = f.(x)
+# Plot the results and compare to target distribution
+x = range(-10.0, 10.0, length=100)
+envelop = [eval_envelop(sampler.envelop, xi) for xi in x]
+target = [f(xi) for xi in x]
 
 histogram(sim, normalize = true, label = "Histogram")
 plot!(x, [target envelop], width = 2, label = ["Normal(μ, σ)" "Envelop"])
@@ -56,12 +57,14 @@ plot!(x, [target envelop], width = 2, label = ["Normal(μ, σ)" "Envelop"])
 f(x) = β^α * x^(α-1) * exp(-β*x) / gamma(α)
 support = (0.0, Inf)
 
+# Build the sampler and simulate 10,000 samples
 sampler = RejectionSampler(f, support)
 @time sim = run_sampler!(sampler, 10000) 
 
-x = linspace(0.0, 5.0, 100)
-envelop = eval_envelop.(sampler.envelop, x)
-target = f.(x)
+# Plot the results and compare to target distribution
+x = range(0.0, 5.0, length=100)
+envelop = [eval_envelop(sampler.envelop, xi) for xi in x]
+target = [f(xi) for xi in x]
 
 histogram(sim, normalize = true, label = "Histogram")
 plot!(x, [target envelop], width = 2, label = ["Gamma(α, β)" "Envelop"])
@@ -71,8 +74,6 @@ plot!(x, [target envelop], width = 2, label = ["Gamma(α, β)" "Envelop"])
 
 
 
-
-
 ![](img/example2.png)
 
 ## Truncated distributions and unknown normalisation constant
@@ -85,12 +86,14 @@ We don't to provide an exact density--it will sample up to proportionality--and
 f(x) = β^α * x^(α-1) * exp(-β*x) / gamma(α)
 support = (1.0, 3.5)
 
+# Build the sampler and simulate 10,000 samples
 sampler = RejectionSampler(f, support)
 @time sim = run_sampler!(sampler, 10000) 
 
-x = linspace(0.01, 8.0, 100)
-envelop = eval_envelop.(sampler.envelop, x)
-target = f.(x)
+# Plot the results and compare to target distribution
+x = range(0.01, 8.0, length=100)
+envelop = [eval_envelop(sampler.envelop, xi) for xi in x]
+target = [f(xi) for xi in x]
 
 histogram(sim, normalize = true, label = "histogram")
 plot!(x, [target envelop], width = 2, label = ["target density" "envelop"])
@@ -99,7 +102,33 @@ plot!(x, [target envelop], width = 2, label = ["target density" "envelop"])
       0.007766 seconds (181.82 k allocations: 3.024 MiB)
 
 
+![](img/example3.png)
+
+## Elastic Net distribution
+
+The following example arises from elastic net regression and smoothing problems. In these cases, the integration constants are not available analytically.
+
+```julia
+# Define function to be sampled
+function f(x, μ, λ1, λ2)
+      δ = x - μ
+      nl = λ1 * abs(δ) + λ2 * δ^2
+      return exp(-nl)
+end
+support = (-Inf, Inf)
+
+# Build the sampler and simulate 10,000 samples
+μ, λ1, λ2 = 0.0, 2.0, 1.0
+sampler = RejectionSampler(x -> f(x, μ, λ1, λ2), support, max_segments = 5)
+@time sim = run_sampler!(sampler, 10000);
 
+# Plot the results and compare to target distribution
+x = range(-2.3, 2.3, length=100)
+envelop = [eval_envelop(sampler.envelop, xi) for xi in x]
+target = [f(xi, μ, λ1, λ2) for xi in x]
 
+histogram(sim, normalize = true, label = "histogram")
+plot!(x, [target envelop], width = 2, label = ["target density" "envelop"])
+```
 
-![](img/example3.png)
+![](img/example4.png)

REQUIRE

@@ -1,3 +1,3 @@
-julia 0.6.0
-ForwardDiff 0.5.0
-StatsBase 0.19.2
+julia 0.7
+ForwardDiff 0.10.1
+StatsBase 0.26.0

src/AdaptiveRejectionSampling.jl

@@ -27,13 +27,13 @@ Finds the horizontal coordinate of the intersection between lines
 """
 function intersection(l1::Line, l2::Line)
     @assert l1.slope != l2.slope "slopes should be different"
-     - (l2.intercept - l1.intercept) / (l2.slope - l1.slope)
+    - (l2.intercept - l1.intercept) / (l2.slope - l1.slope)
 end
 
 """
     exp_integral(l::Line, x1::Float64, x2::Float64)
 Computes the integral
-    ``LaTeX \int_{x_1} ^ {x_2} \exp\{ax + b\} dx. ``
+    ``LaTeX \\int_{x_1} ^ {x_2} \\exp\\{ax + b\\} dx. ``
 The resulting value is the weight assigned to the segment [x1, x2] in the envelop
 """
 function exp_integral(l::Line, x1::Float64, x2::Float64)
@@ -50,8 +50,8 @@ over c_1 to c_k+1 is one, so that the envelop is interpreted as a density.
 """
 mutable struct Envelop
     lines::Vector{Line}
-    cutpoints::Vector{Float64}
-    weights::Vector{Float64}
+    cutpoints::AbstractVector{Float64}
+    weights::AbstractVector{Float64}
     size::Int
 
     Envelop(lines::Vector{Line}, support::Tuple{Float64, Float64}) = begin
@@ -161,7 +161,7 @@ interval in which f has positive value, and zero elsewhere.
 - `max_failed_factor::Float64 = 0.001`: level at which throw an error if one single sample has a rejection rate
     exceeding this value
 """
-mutable struct RejectionSampler
+struct RejectionSampler
     objective::Objective
     envelop::Envelop
     max_segments::Int
@@ -201,7 +201,7 @@ mutable struct RejectionSampler
         grid_lims = max(search_range[1], support[1]), min(search_range[2], support[2])
         grid = grid_lims[1]:δ:grid_lims[2]
         i1, i2 = findfirst(grad.(grid) .> 0.), findfirst(grad.(grid) .< 0.)
-        @assert (i1 > 0) &&  (i2 > 0) "couldn't find initial points, please provide them or change `search_range`"
+        @assert (i1 != nothing) &&  (i2 != nothing) "couldn't find initial points, please provide them or change `search_range`"
         x1, x2 = grid[i1], grid[i2]
         RejectionSampler(f, support, (x1, x2); kwargs...)
     end

test/runtests.jl

@@ -1,7 +1,7 @@
 import AdaptiveRejectionSampling
 ars = AdaptiveRejectionSampling
 
-using Base.Test
+using Test
 
 @testset "Line" begin
     @test ars.Line(2.0, 3) isa ars.Line