From 43dabf5626a98900874041c60b864009f74f8a6a Mon Sep 17 00:00:00 2001 From: Tor Erlend Fjelde Date: Mon, 20 Dec 2021 00:22:51 +0000 Subject: [PATCH 01/13] pass model around instead of prior_transform and loglike --- src/NestedSamplers.jl | 13 +++-- src/model.jl | 38 ++++++++++++- src/proposals/Proposals.jl | 114 +++++++++++++++++++------------------ src/step.jl | 22 +++---- 4 files changed, 112 insertions(+), 75 deletions(-) diff --git a/src/NestedSamplers.jl b/src/NestedSamplers.jl index 42262ed8..1c631fa8 100644 --- a/src/NestedSamplers.jl +++ b/src/NestedSamplers.jl @@ -1,11 +1,5 @@ module NestedSamplers -# load submodules -include("bounds/Bounds.jl") -using .Bounds -include("proposals/Proposals.jl") -using .Proposals - using LinearAlgebra using Random using Random: AbstractRNG, GLOBAL_RNG @@ -32,6 +26,13 @@ export Bounds, Nested include("model.jl") # The default model for nested sampling + +# load submodules +include("bounds/Bounds.jl") +using .Bounds +include("proposals/Proposals.jl") +using .Proposals + include("staticsampler.jl") # The static nested sampler include("step.jl") # The stepping mechanics (extends AbstractMCMC) include("sample.jl") # Custom sampling (extends AbstractMCMC) diff --git a/src/model.jl b/src/model.jl index 5de844f8..409b6aec 100644 --- a/src/model.jl +++ b/src/model.jl @@ -1,3 +1,17 @@ +struct PriorTransformAndLogLike{T,L} + prior_transform::T + loglikelihood::L +end + +function (f::PriorTransformAndLogLike)(args...) + return f.prior_transform(args...), f.loglikelihood(args...) +end + +prior_transform(f::PriorTransformAndLogLike, args...) = f.prior_transform(args...) +loglikelihood(f::PriorTransformAndLogLike, args...) = f.loglikelihood(args...) +function prior_transform_and_loglikelihood(f::PriorTransformAndLogLike, args...) + return f.prior_transform(args...), f.loglikelihood(args...) +end """ NestedModel(loglike, prior_transform) @@ -10,12 +24,30 @@ **Note:** `loglike` is the only function used for likelihood calculations. This means if you want your priors to be used for the likelihood calculations they must be manually included in the `loglike` function. """ -struct NestedModel <: AbstractModel - loglike::Function - prior_transform::Function +struct NestedModel{F} <: AbstractModel + prior_transform_and_loglike::F +end + +function NestedModel(loglike, prior_transform) + return NestedModel(PriorTransformAndLogLike(prior_transform, loglike)) end function NestedModel(loglike, priors::AbstractVector{<:UnivariateDistribution}) prior_transform(X) = quantile.(priors, X) return NestedModel(loglike, prior_transform) end + +function prior_transform(model::NestedModel{<:PriorTransformAndLogLike}, args...) + return prior_transform(model.prior_transform_and_loglike, args...) +end + +function loglikelihood(model::NestedModel{<:PriorTransformAndLogLike}, args...) + return loglikelihood(model.prior_transform_and_loglike, args...) +end + +function prior_transform_and_loglikelihood( + model::NestedModel{<:PriorTransformAndLogLike}, + args... +) + return prior_transform_and_loglikelihood(model.prior_transform_and_loglike, args...) +end diff --git a/src/proposals/Proposals.jl b/src/proposals/Proposals.jl index 4c83657c..1fc429ca 100644 --- a/src/proposals/Proposals.jl +++ b/src/proposals/Proposals.jl @@ -12,6 +12,7 @@ The available implementations are """ module Proposals +using ..NestedSamplers: prior_transform_and_loglikelihood using ..Bounds using Random @@ -54,19 +55,19 @@ end @deprecate Uniform() Rejection() -function (prop::Rejection)(rng::AbstractRNG, +function (prop::Rejection)( + rng::AbstractRNG, point::AbstractVector, logl_star, bounds::AbstractBoundingSpace, - loglike, - prior_transform) + model +) ncall = 0 for _ in 1:prop.maxiter u = rand(rng, bounds) unitcheck(u) || continue - v = prior_transform(u) - logl = loglike(v) + v, logl = prior_transform_and_loglikelihood(model, u) ncall += 1 logl ≥ logl_star && return u, v, logl, ncall end @@ -95,13 +96,14 @@ Propose a new live point by random walking away from an existing live point. @assert scale ≥ 0 "Proposal scale must be non-negative" end -function (prop::RWalk)(rng::AbstractRNG, - point::AbstractVector, - logl_star, - bounds::AbstractBoundingSpace, - loglike, - prior_transform; - kwargs...) +function (prop::RWalk)( + rng::AbstractRNG, + point::AbstractVector, + logl_star, + bounds::AbstractBoundingSpace, + model; + kwargs... +) # setup n = length(point) scale_init = prop.scale @@ -129,8 +131,7 @@ function (prop::RWalk)(rng::AbstractRNG, end end # check proposed point - v_prop = prior_transform(u_prop) - logl_prop = loglike(v_prop) + v_prop, logl_prop = prior_transfrom_and_loglike(model, u_prop) if logl_prop ≥ logl_star u = u_prop v = v_prop @@ -188,13 +189,14 @@ proposals. @assert scale ≥ 0 "Proposal scale must be non-negative" end -function (prop::RStagger)(rng::AbstractRNG, - point::AbstractVector, - logl_star, - bounds::AbstractBoundingSpace, - loglike, - prior_transform; - kwargs...) +function (prop::RStagger)( + rng::AbstractRNG, + point::AbstractVector, + logl_star, + bounds::AbstractBoundingSpace, + model; + kwargs... +) #setup n = length(point) scale_init = prop.scale @@ -223,8 +225,7 @@ function (prop::RStagger)(rng::AbstractRNG, end end # check proposed point - v_prop = prior_transform(u_prop) - logl_prop = loglike(v_prop) + v_prop, logl_prop = prior_transform_and_loglikelihood(model, u_prop) if logl_prop ≥ logl_star u = u_prop v = v_prop @@ -276,13 +277,14 @@ This is a standard _Gibbs-like_ implementation where a single multivariate slice @assert scale ≥ 0 "Proposal scale must be non-negative" end -function (prop::Slice)(rng::AbstractRNG, - point::AbstractVector, - logl_star, - bounds::AbstractBoundingSpace, - loglike, - prior_transform; - kwargs...) +function (prop::Slice)( + rng::AbstractRNG, + point::AbstractVector, + logl_star, + bounds::AbstractBoundingSpace, + model; + kwargs... +) # setup n = length(point) nc = nexpand = ncontract = 0 @@ -303,8 +305,11 @@ function (prop::Slice)(rng::AbstractRNG, # select axis axis = axes[idx, :] - u, v, logl, nc, nexpand, ncontract = sample_slice(rng, axis, point, logl_star, loglike, - prior_transform, nc, nexpand, ncontract) + u, v, logl, nc, nexpand, ncontract = sample_slice( + rng, axis, point, logl_star, + model, + nc, nexpand, ncontract + ) end # end of slice sample along a random direction end # end of slice sampling loop @@ -330,13 +335,14 @@ This is a standard _random_ implementation where each slice is along a random di @assert scale ≥ 0 "Proposal scale must be non-negative" end -function (prop::RSlice)(rng::AbstractRNG, - point::AbstractVector, - logl_star, - bounds::AbstractBoundingSpace, - loglike, - prior_transform; - kwargs...) +function (prop::RSlice)( + rng::AbstractRNG, + point::AbstractVector, + logl_star, + bounds::AbstractBoundingSpace, + model; + kwargs... +) # setup n = length(point) nc = nexpand = ncontract = 0 @@ -350,8 +356,11 @@ function (prop::RSlice)(rng::AbstractRNG, # transform and scale into parameter space axis = prop.scale .* (Bounds.axes(bounds) * drhat) - u, v, logl, nc, nexpand, ncontract = sample_slice(rng, axis, point, logl_star, loglike, - prior_transform, nc, nexpand, ncontract) + u, v, logl, nc, nexpand, ncontract = sample_slice( + rng, axis, point, logl_star, + model, + nc, nexpand, ncontract + ) end # end of random slice sampling loop # update random slice proposal scale based on the relative size of the slices compared to the initial guess @@ -361,13 +370,12 @@ function (prop::RSlice)(rng::AbstractRNG, end # end of function RSlice # Method for slice sampling -function sample_slice(rng, axis, u, logl_star, loglike, prior_transform, nc, nexpand, ncontract) +function sample_slice(rng, axis, u, logl_star, model, nc, nexpand, ncontract) # define starting window r = rand(rng) # initial scale/offset u_l = @. u - r * axis # left bound if unitcheck(u_l) - v_l = prior_transform(u_l) - logl_l = loglike(v_l) + v_l, logl_l = prior_transform_and_loglikelihood(model, u_l) else logl_l = -Inf end @@ -376,8 +384,7 @@ function sample_slice(rng, axis, u, logl_star, loglike, prior_transform, nc, nex u_r = u_l .+ axis # right bound if unitcheck(u_r) - v_r = prior_transform(u_r) - logl_r = loglike(v_r) + v_r, logl_r = prior_transform_and_loglikelihood(model, u_r) else logl_r = -Inf end @@ -387,9 +394,8 @@ function sample_slice(rng, axis, u, logl_star, loglike, prior_transform, nc, nex # stepping out left and right bounds while logl_l ≥ logl_star u_l .-= axis - if unitcheck(u_l) - v_l = prior_transform(u_l) - logl_l = loglike(v_l) + if unitcheck(u_l) + v_l, logl_l = prior_transform_and_loglikelihood(model, u_l) else logl_l = -Inf end @@ -399,9 +405,8 @@ function sample_slice(rng, axis, u, logl_star, loglike, prior_transform, nc, nex while logl_r ≥ logl_star u_r .+= axis - if unitcheck(u_r) - v_r = prior_transform(u_r) - logl_r = loglike(v_r) + if unitcheck(u_r) + v_r, logl_r = prior_transform_and_loglikelihood(model, u_r) else logl_r = -Inf end @@ -422,9 +427,8 @@ function sample_slice(rng, axis, u, logl_star, loglike, prior_transform, nc, nex # propose a new position r = rand(rng) u_prop = @. u_l + r * u_hat - if unitcheck(u_prop) - v_prop = prior_transform(u_prop) - logl_prop = loglike(v_prop) + if unitcheck(u_prop) + u_prop, logl_prop = prior_transform_and_loglikelihood(model, u_prop) else logl_prop = -Inf end diff --git a/src/step.jl b/src/step.jl index 62729b49..55b5a438 100644 --- a/src/step.jl +++ b/src/step.jl @@ -17,7 +17,7 @@ function step(rng, model, sampler::Nested; kwargs...) point = rand(rng, eltype(us), sampler.ndims) bound = Bounds.fit(Bounds.NoBounds, us) proposal = Proposals.Rejection() - u, v, ll, nc = proposal(rng, v_dead, logl_dead, bound, model.loglike, model.prior_transform) + u, v, ll, nc = proposal(rng, v_dead, logl_dead, bound, model) us[:, idx_dead] .= u vs[:, idx_dead] .= v @@ -75,12 +75,12 @@ function step(rng, model, sampler, state; kwargs...) since_update = 0 point, bound = rand_live(rng, active_bound, point[:, :]) end - u, v, logl, nc = sampler.proposal(rng, point, logl_dead, bound, model.loglike, model.prior_transform) + u, v, logl, nc = sampler.proposal(rng, point, logl_dead, bound, model) else point = rand(rng, eltype(state.us), sampler.ndims) bound = Bounds.fit(Bounds.NoBounds, state.us) proposal = Proposals.Rejection() - u, v, logl, nc = proposal(rng, point, logl_dead, bound, model.loglike, model.prior_transform) + u, v, logl, nc = proposal(rng, point, logl_dead, bound, model) end state.us[:, idx_dead] .= u @@ -167,24 +167,24 @@ end ## Helpers -init_particles(rng, ndims, nactive, prior, loglike) = - init_particles(rng, Float64, ndims, nactive, prior, loglike) +init_particles(rng, ndims, nactive, model) = + init_particles(rng, Float64, ndims, nactive, model) init_particles(rng, model, sampler) = - init_particles(rng, sampler.ndims, sampler.nactive, model.prior_transform, model.loglike) + init_particles(rng, sampler.ndims, sampler.nactive, model) # loop and fill arrays, checking validity of points # will retry 100 times before erroring -function init_particles(rng, T, ndims, nactive, prior, loglike) +function init_particles(rng, T, ndims, nactive, model) us = rand(rng, T, ndims, nactive) - vs = mapslices(prior, us, dims=1) - logl = dropdims(mapslices(loglike, vs, dims=1), dims=1) + vs = mapslices(Base.Fix1(prior_transform, model), us, dims=1) + logl = dropdims(mapslices(Base.Fix1(loglikelihood, model), vs, dims=1), dims=1) ntries = 1 while true any(isfinite, logl) && break us .= rand(rng, T, ndims, nactive) - vs .= mapslices(prior, us, dims=1) - logl .= mapslices(loglike, vs, dims=1) + vs .= mapslices(Base.Fix1(prior_transform, model), us, dims=1) + logl .= mapslices(Base.Fix1(loglikelihood, model), vs, dims=1) ntries += 1 ntries > 100 && error("After 100 attempts, could not initialize any live points with finite loglikelihood. Please check your prior transform and loglikelihood methods.") end From c7a7c38ee4b08aa2dfd8e16af618afb5590cfae9 Mon Sep 17 00:00:00 2001 From: Tor Erlend Fjelde Date: Mon, 20 Dec 2021 00:23:49 +0000 Subject: [PATCH 02/13] update README --- README.md | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/README.md b/README.md index 12f114bc..06b2c420 100644 --- a/README.md +++ b/README.md @@ -29,18 +29,19 @@ For in-depth usage, see the [online documentation](https://turinglang.github.io/ ```julia using NestedSamplers using Distributions +using LinearAlgebra -logl(X) = exp(-(X - [1, -1]) / 2) +logl(X) = logpdf(MvNormal([1, -1], I), X) prior(X) = 4 .* (X .- 0.5) # or equivalently -prior = [Uniform(-2, 2), Uniform(-2, 2)] -model = NestedModel(logl, prior) +priors = [Uniform(-2, 2), Uniform(-2, 2)] +model = NestedModel(logl, priors) ``` after defining the model, set up the nested sampler. This will involve choosing the bounding space and proposal scheme, or you can rely on the defaults. In addition, we need to define the dimensionality of the problem and the number of live points. More points results in a more precise evidence estimate at the cost of runtime. For more information, see the docs. ```julia -bounds = Bounds.MultiElliipsoid +bounds = Bounds.MultiEllipsoid prop = Proposals.Slice(slices=10) # 1000 live points sampler = Nested(2, 1000; bounds=bounds, proposal=prop) From 61e97c39fd2972b2b89506edef15cf529821dff2 Mon Sep 17 00:00:00 2001 From: Tor Erlend Fjelde Date: Mon, 20 Dec 2021 00:29:12 +0000 Subject: [PATCH 03/13] added some default implementations --- src/model.jl | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) diff --git a/src/model.jl b/src/model.jl index 409b6aec..4c2887bc 100644 --- a/src/model.jl +++ b/src/model.jl @@ -37,17 +37,22 @@ function NestedModel(loglike, priors::AbstractVector{<:UnivariateDistribution}) return NestedModel(loglike, prior_transform) end +function prior_transform(model, args...) + return first(prior_transform_and_loglikelihood(model, args...)) +end + function prior_transform(model::NestedModel{<:PriorTransformAndLogLike}, args...) return prior_transform(model.prior_transform_and_loglike, args...) end +function loglikelihood(model, args...) + return last(prior_transform_and_loglikelihood(model, args...)) +end + function loglikelihood(model::NestedModel{<:PriorTransformAndLogLike}, args...) return loglikelihood(model.prior_transform_and_loglike, args...) end -function prior_transform_and_loglikelihood( - model::NestedModel{<:PriorTransformAndLogLike}, - args... -) +function prior_transform_and_loglikelihood(model::NestedModel, args...) return prior_transform_and_loglikelihood(model.prior_transform_and_loglike, args...) end From ef1f0da2e5b561cb528293e3e8f5336e44d91eca Mon Sep 17 00:00:00 2001 From: Tor Erlend Fjelde Date: Mon, 20 Dec 2021 00:53:49 +0000 Subject: [PATCH 04/13] assume by default that prior_transform_and_loglike is a function --- src/model.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/model.jl b/src/model.jl index 4c2887bc..4327f650 100644 --- a/src/model.jl +++ b/src/model.jl @@ -54,5 +54,5 @@ function loglikelihood(model::NestedModel{<:PriorTransformAndLogLike}, args...) end function prior_transform_and_loglikelihood(model::NestedModel, args...) - return prior_transform_and_loglikelihood(model.prior_transform_and_loglike, args...) + return model.prior_transform_and_loglike(args...) end From c9246f6a600d34a3044c62ee6de612e415952ef4 Mon Sep 17 00:00:00 2001 From: Tor Erlend Fjelde Date: Mon, 20 Dec 2021 01:06:07 +0000 Subject: [PATCH 05/13] fixed init_particles --- src/step.jl | 19 ++++++++++++++----- 1 file changed, 14 insertions(+), 5 deletions(-) diff --git a/src/step.jl b/src/step.jl index 55b5a438..24725c52 100644 --- a/src/step.jl +++ b/src/step.jl @@ -177,14 +177,23 @@ init_particles(rng, model, sampler) = # will retry 100 times before erroring function init_particles(rng, T, ndims, nactive, model) us = rand(rng, T, ndims, nactive) - vs = mapslices(Base.Fix1(prior_transform, model), us, dims=1) - logl = dropdims(mapslices(Base.Fix1(loglikelihood, model), vs, dims=1), dims=1) - ntries = 1 + vs_and_logl = mapslices( + Base.Fix1(prior_transform_and_loglikelihood, model), us; + dims=1 + ) + vs = map(first, vs_and_logl) + logl = map(last, vs_and_logl) + + ntries = 1 while true any(isfinite, logl) && break us .= rand(rng, T, ndims, nactive) - vs .= mapslices(Base.Fix1(prior_transform, model), us, dims=1) - logl .= mapslices(Base.Fix1(loglikelihood, model), vs, dims=1) + vs_and_logl .= mapslices( + Base.Fix1(prior_transform_and_loglikelihood, model), us; + dims=1 + ) + vs. = map(first, vs_and_logl) + logl. = map(last, vs_and_logl) ntries += 1 ntries > 100 && error("After 100 attempts, could not initialize any live points with finite loglikelihood. Please check your prior transform and loglikelihood methods.") end From 0ebff71c7c3152afbfc2ebf7f5b58ffea14273c0 Mon Sep 17 00:00:00 2001 From: Tor Erlend Fjelde Date: Mon, 20 Dec 2021 01:39:08 +0000 Subject: [PATCH 06/13] fixed init_particles for real --- src/step.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/step.jl b/src/step.jl index 24725c52..8041bb17 100644 --- a/src/step.jl +++ b/src/step.jl @@ -181,19 +181,19 @@ function init_particles(rng, T, ndims, nactive, model) Base.Fix1(prior_transform_and_loglikelihood, model), us; dims=1 ) - vs = map(first, vs_and_logl) - logl = map(last, vs_and_logl) + vs = mapreduce(first, hcat, vs_and_logl) + logl = dropdims(map(last, vs_and_logl), dims=1) ntries = 1 while true any(isfinite, logl) && break - us .= rand(rng, T, ndims, nactive) + rand!(rng, us) vs_and_logl .= mapslices( Base.Fix1(prior_transform_and_loglikelihood, model), us; dims=1 ) - vs. = map(first, vs_and_logl) - logl. = map(last, vs_and_logl) + vs .= mapreduce(first, hcat, vs_and_logl) + map!(last, logl, vs_and_logl) ntries += 1 ntries > 100 && error("After 100 attempts, could not initialize any live points with finite loglikelihood. Please check your prior transform and loglikelihood methods.") end From b65a6603bc8eb51ea39f0c2c318663d17d42a4bd Mon Sep 17 00:00:00 2001 From: Tor Erlend Fjelde Date: Mon, 20 Dec 2021 15:56:13 +0000 Subject: [PATCH 07/13] bump minor version --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index add615ec..3d0e2262 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "NestedSamplers" uuid = "41ceaf6f-1696-4a54-9b49-2e7a9ec3782e" authors = ["Miles Lucas "] -version = "0.7.1" +version = "0.8.0" [deps] AbstractMCMC = "80f14c24-f653-4e6a-9b94-39d6b0f70001" From bcf310d1905846a57ea62bf2fcb590cc242da430 Mon Sep 17 00:00:00 2001 From: Tor Erlend Fjelde Date: Tue, 21 Dec 2021 01:21:09 +0000 Subject: [PATCH 08/13] fix tests --- test/proposals/proposals.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/proposals/proposals.jl b/test/proposals/proposals.jl index f3e72800..14529da1 100644 --- a/test/proposals/proposals.jl +++ b/test/proposals/proposals.jl @@ -18,7 +18,7 @@ const BOUNDS = [ us = 0.7 .* rand(rng, 2, 10) # live points should be within the ellipsoid point, _bound = Bounds.rand_live(rng, bound, us) loglstar = logl(prior(point)) - u, v, logL = prop(rng, point, loglstar, _bound, logl, prior) + u, v, logL = prop(rng, point, loglstar, _bound, NestedSamplers.PriorTransformAndLogLikelihood(prior, logl)) # simple bounds checks @test all(x -> 0 < x < 1, u) @test all(x -> -1 < x < 1, v) @@ -82,4 +82,4 @@ end @test_throws AssertionError Proposals.RSlice(slices=-2) @test_throws AssertionError Proposals.RSlice(scale=-3) -end \ No newline at end of file +end From 29ed709ee289b1893be66b949e93c20b67196f54 Mon Sep 17 00:00:00 2001 From: Tor Erlend Fjelde Date: Tue, 21 Dec 2021 01:25:37 +0000 Subject: [PATCH 09/13] rename PriorTransformAndLogLike to PriorTransformAndLogLikelihood --- src/model.jl | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/model.jl b/src/model.jl index 4327f650..1ad841e2 100644 --- a/src/model.jl +++ b/src/model.jl @@ -1,15 +1,15 @@ -struct PriorTransformAndLogLike{T,L} +struct PriorTransformAndLogLikelihood{T,L} prior_transform::T loglikelihood::L end -function (f::PriorTransformAndLogLike)(args...) +function (f::PriorTransformAndLogLikelihood)(args...) return f.prior_transform(args...), f.loglikelihood(args...) end -prior_transform(f::PriorTransformAndLogLike, args...) = f.prior_transform(args...) -loglikelihood(f::PriorTransformAndLogLike, args...) = f.loglikelihood(args...) -function prior_transform_and_loglikelihood(f::PriorTransformAndLogLike, args...) +prior_transform(f::PriorTransformAndLogLikelihood, args...) = f.prior_transform(args...) +loglikelihood(f::PriorTransformAndLogLikelihood, args...) = f.loglikelihood(args...) +function prior_transform_and_loglikelihood(f::PriorTransformAndLogLikelihood, args...) return f.prior_transform(args...), f.loglikelihood(args...) end @@ -29,7 +29,7 @@ struct NestedModel{F} <: AbstractModel end function NestedModel(loglike, prior_transform) - return NestedModel(PriorTransformAndLogLike(prior_transform, loglike)) + return NestedModel(PriorTransformAndLogLikelihood(prior_transform, loglike)) end function NestedModel(loglike, priors::AbstractVector{<:UnivariateDistribution}) @@ -41,7 +41,7 @@ function prior_transform(model, args...) return first(prior_transform_and_loglikelihood(model, args...)) end -function prior_transform(model::NestedModel{<:PriorTransformAndLogLike}, args...) +function prior_transform(model::NestedModel{<:PriorTransformAndLogLikelihood}, args...) return prior_transform(model.prior_transform_and_loglike, args...) end @@ -49,7 +49,7 @@ function loglikelihood(model, args...) return last(prior_transform_and_loglikelihood(model, args...)) end -function loglikelihood(model::NestedModel{<:PriorTransformAndLogLike}, args...) +function loglikelihood(model::NestedModel{<:PriorTransformAndLogLikelihood}, args...) return loglikelihood(model.prior_transform_and_loglike, args...) end From 644fb50417fab1c2c5c617559d02b1fc1e7f76d9 Mon Sep 17 00:00:00 2001 From: Tor Erlend Fjelde Date: Tue, 21 Dec 2021 01:38:55 +0000 Subject: [PATCH 10/13] fix for PriorTransformAndLogLikelihood --- src/model.jl | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/src/model.jl b/src/model.jl index 1ad841e2..bf1b46f7 100644 --- a/src/model.jl +++ b/src/model.jl @@ -3,15 +3,16 @@ struct PriorTransformAndLogLikelihood{T,L} loglikelihood::L end -function (f::PriorTransformAndLogLikelihood)(args...) - return f.prior_transform(args...), f.loglikelihood(args...) +function (f::PriorTransformAndLogLikelihood)(u) + v = f.prior_transform(u) + return (v, f.loglikelihood(v)) end -prior_transform(f::PriorTransformAndLogLikelihood, args...) = f.prior_transform(args...) -loglikelihood(f::PriorTransformAndLogLikelihood, args...) = f.loglikelihood(args...) -function prior_transform_and_loglikelihood(f::PriorTransformAndLogLikelihood, args...) - return f.prior_transform(args...), f.loglikelihood(args...) +prior_transform(f::PriorTransformAndLogLikelihood, u) = f.prior_transform(u) +function loglikelihood_from_uniform(f::PriorTransformAndLogLikelihood, u) + return last(prior_transform_and_loglikelihood(f, u)) end +prior_transform_and_loglikelihood(f::PriorTransformAndLogLikelihood, u) = f(u) """ NestedModel(loglike, prior_transform) From e8bc9ee21740f6353c91038435f74c853b3ea125 Mon Sep 17 00:00:00 2001 From: Tor Erlend Fjelde Date: Tue, 21 Dec 2021 01:39:21 +0000 Subject: [PATCH 11/13] renamed loglikelihood method to loglikelihood_from_uniform --- src/model.jl | 12 ++++++------ src/proposals/Proposals.jl | 2 +- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/src/model.jl b/src/model.jl index bf1b46f7..5ab9d446 100644 --- a/src/model.jl +++ b/src/model.jl @@ -26,7 +26,7 @@ prior_transform_and_loglikelihood(f::PriorTransformAndLogLikelihood, u) = f(u) `loglike` is the only function used for likelihood calculations. This means if you want your priors to be used for the likelihood calculations they must be manually included in the `loglike` function. """ struct NestedModel{F} <: AbstractModel - prior_transform_and_loglike::F + prior_transform_and_loglikelihood::F end function NestedModel(loglike, prior_transform) @@ -43,17 +43,17 @@ function prior_transform(model, args...) end function prior_transform(model::NestedModel{<:PriorTransformAndLogLikelihood}, args...) - return prior_transform(model.prior_transform_and_loglike, args...) + return prior_transform(model.prior_transform_and_loglikelihood, args...) end -function loglikelihood(model, args...) +function loglikelihood_from_uniform(model, args...) return last(prior_transform_and_loglikelihood(model, args...)) end -function loglikelihood(model::NestedModel{<:PriorTransformAndLogLikelihood}, args...) - return loglikelihood(model.prior_transform_and_loglike, args...) +function loglikelihood_from_uniform(model::NestedModel{<:PriorTransformAndLogLikelihood}, args...) + return loglikelihood_from_uniform(model.prior_transform_and_loglikelihood, args...) end function prior_transform_and_loglikelihood(model::NestedModel, args...) - return model.prior_transform_and_loglike(args...) + return model.prior_transform_and_loglikelihood(args...) end diff --git a/src/proposals/Proposals.jl b/src/proposals/Proposals.jl index 1fc429ca..021fdbd6 100644 --- a/src/proposals/Proposals.jl +++ b/src/proposals/Proposals.jl @@ -131,7 +131,7 @@ function (prop::RWalk)( end end # check proposed point - v_prop, logl_prop = prior_transfrom_and_loglike(model, u_prop) + v_prop, logl_prop = prior_transform_and_loglikelihood(model, u_prop) if logl_prop ≥ logl_star u = u_prop v = v_prop From 2b129c0325c625a085f9871ae6bf41de63fcfefb Mon Sep 17 00:00:00 2001 From: Tor Erlend Fjelde Date: Tue, 21 Dec 2021 03:03:14 +0000 Subject: [PATCH 12/13] fixed a bug with Slice proposal --- src/proposals/Proposals.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/proposals/Proposals.jl b/src/proposals/Proposals.jl index 021fdbd6..42954bd6 100644 --- a/src/proposals/Proposals.jl +++ b/src/proposals/Proposals.jl @@ -428,7 +428,7 @@ function sample_slice(rng, axis, u, logl_star, model, nc, nexpand, ncontract) r = rand(rng) u_prop = @. u_l + r * u_hat if unitcheck(u_prop) - u_prop, logl_prop = prior_transform_and_loglikelihood(model, u_prop) + v_prop, logl_prop = prior_transform_and_loglikelihood(model, u_prop) else logl_prop = -Inf end From bfb05c8e8ba770d955288d6c993eb6e96607c84b Mon Sep 17 00:00:00 2001 From: Tor Erlend Fjelde Date: Tue, 21 Dec 2021 21:27:16 +0000 Subject: [PATCH 13/13] add rng to lines in tests which are missing rng --- test/models.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/models.jl b/test/models.jl index 2a78cfa1..7d97b072 100644 --- a/test/models.jl +++ b/test/models.jl @@ -9,7 +9,7 @@ const test_props = [Proposals.Rejection(maxiter=Int(1e6)), Proposals.RWalk(ratio sampler = Nested(D, 50D; bounds=bound, proposal=proposal) chain, state = sample(rng, model, sampler; dlogz=0.01) - chain_res = sample(chain, Weights(vec(chain[:weights])), length(chain)) + chain_res = sample(rng, chain, Weights(vec(chain[:weights])), length(chain)) # test posteriors vals = Array(chain_res) means = mean(vals, dims=1) @@ -64,7 +64,7 @@ const test_props = [Proposals.Rejection(maxiter=Int(1e6)), Proposals.RWalk(ratio spl = Nested(2, 1000, bounds=bound, proposal=proposal) chain, state = sample(rng, model, spl; dlogz=0.01) - chain_res = sample(chain, Weights(vec(chain[:weights])), length(chain)) + chain_res = sample(rng, chain, Weights(vec(chain[:weights])), length(chain)) diff = state.logz - analytic_logz atol = 6state.logzerr @@ -90,7 +90,7 @@ const test_props = [Proposals.Rejection(maxiter=Int(1e6)), Proposals.RWalk(ratio @test state.logz ≈ logz atol = 5state.logzerr - chain_res = sample(chain, Weights(vec(chain[:weights])), length(chain)) + chain_res = sample(rng, chain, Weights(vec(chain[:weights])), length(chain)) xmodes = sort!(findpeaks(chain_res[:, 1, 1])[1:5]) @test all(isapprox.(xmodes, 0.1:0.2:0.9, atol=0.2)) ymodes = sort!(findpeaks(chain_res[:, 2, 1])[1:5])