Merge pull request #30 from devmotion/fix

Fix emcee

Author: cpfiffer

src/emcee.jl

@@ -90,7 +90,7 @@ function move(
     alphamult = (n - 1) * log(z)
 
     # Make new parameters
-    y = @. walker.params + z * (other_walker.params - walker.params)
+    y = @. other_walker.params + z * (walker.params - other_walker.params)
 
     # Construct a new walker
     new_walker = Transition(model, y)
@@ -215,4 +215,4 @@ end
 #             theta = theta_min + (theta_max - theta_min) * rand()
 #         end
 #     end 
-# end
+# end

test/emcee.jl

@@ -0,0 +1,52 @@
+@testset "emcee.jl" begin
+    @testset "example" begin
+        @testset "untransformed space" begin
+            # define model
+            function logprob(θ)
+                s, m = θ
+                s > 0 || return -Inf
+
+                mdist = Normal(0, sqrt(s))
+                obsdist = Normal(m, sqrt(s))
+
+                return logpdf(InverseGamma(2, 3), s) + logpdf(mdist, m) +
+                    logpdf(obsdist, 1.5) + logpdf(obsdist, 2.0)
+            end
+            model = DensityModel(logprob)
+
+            # perform stretch move and sample from prior in initial step
+            Random.seed!(100)
+            sampler = Ensemble(1_000, StretchProposal([InverseGamma(2, 3), Normal(0, 1)]))
+            chain = sample(model, sampler, 1_000;
+                           param_names = ["s", "m"], chain_type = Chains)
+
+            @test mean(chain["s"].value) ≈ 49/24 atol=0.1
+            @test mean(chain["m"].value) ≈ 7/6 atol=0.1
+        end
+
+        @testset "transformed space" begin
+            # define model
+            function logprob(θ)
+                logs, m = θ
+                s = exp(logs)
+                sqrts = sqrt(s)
+
+                mdist = Normal(0, sqrts)
+                obsdist = Normal(m, sqrts)
+
+                return logpdf(InverseGamma(2, 3), s) + logpdf(mdist, m) +
+                    logpdf(obsdist, 1.5) + logpdf(obsdist, 2.0) + logs
+            end
+            model = DensityModel(logprob)
+
+            # perform stretch move and sample from normal distribution in initial step
+            Random.seed!(100)
+            sampler = Ensemble(1_000, StretchProposal(MvNormal(2, 1)))
+            chain = sample(model, sampler, 1_000;
+                           param_names = ["logs", "m"], chain_type = Chains)
+
+            @test mean(exp.(chain["logs"].value)) ≈ 49/24 atol=0.1
+            @test mean(chain["m"].value) ≈ 7/6 atol=0.1
+        end
+    end
+end

test/runtests.jl

@@ -97,5 +97,7 @@ using Test
 
         @test chain1[1].params == val
     end
+
+    @testset "EMCEE" begin include("emcee.jl") end
 end