-
Notifications
You must be signed in to change notification settings - Fork 17
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Adaptive proposals #39
base: master
Are you sure you want to change the base?
Changes from 21 commits
e9fd602
5f0ddfa
f42b784
a188780
d8989aa
565f12a
cc16195
802ec67
c7623c4
a33937e
4007bd0
71e010b
2d59ede
279aea7
93f17c5
16715e1
046c21b
b91fcc0
387eff4
8fb1048
4071675
a63262d
afd3ed1
fe8562c
f66f647
2988ff2
e5ad041
7aa8631
4999349
d4b3f6b
cb52c7f
5a2a175
b2be967
518aab1
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
There are no files selected for viewing
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,107 @@ | ||
""" | ||
Adaptor(; tune=25, target=0.44, bound=10., δmax=0.2) | ||
|
||
A helper struct for univariate adaptive proposal kernels. This tracks the | ||
number of accepted proposals and the total number of attempted proposals. The | ||
proposal kernel is tuned every `tune` proposals, such that the scale (log(σ) in | ||
the case of a Normal kernel, log(b) for a Uniform kernel) of the proposal is | ||
increased (decreased) by `δ(n) = min(δmax, 1/√n)` at tuning step `n` if the | ||
estimated acceptance probability is higher (lower) than `target`. The target | ||
acceptance probability defaults to 0.44 which is supposedly optimal for 1D | ||
proposals. To ensure ergodicity, the scale of the proposal has to be bounded | ||
(by `bound`), although this is often not required in practice. | ||
""" | ||
mutable struct Adaptor | ||
accepted::Int | ||
total::Int | ||
tune::Int # tuning interval | ||
target::Float64 # target acceptance rate | ||
bound::Float64 # bound on logσ of Gaussian kernel | ||
δmax::Float64 # maximum adaptation step | ||
end | ||
|
||
function Adaptor(; tune=25, target=0.44, bound=10., δmax=0.2) | ||
return Adaptor(0, 0, tune, target, bound, δmax) | ||
end | ||
|
||
""" | ||
AdaptiveProposal{P} | ||
|
||
An adaptive Metropolis-Hastings proposal. In order for this to work, the | ||
proposal kernel should implement the `adapted(proposal, δ)` method, where `δ` | ||
is the increment/decrement applied to the scale of the proposal distribution | ||
during adaptation (e.g. for a Normal distribution the scale is `log(σ)`, so | ||
that after adaptation the proposal is `Normal(0, exp(log(σ) + δ))`). | ||
|
||
# Example | ||
```julia | ||
julia> p = AdaptiveProposal(Uniform(-0.2, 0.2)); | ||
|
||
julia> rand(p) | ||
0.07975590594518434 | ||
``` | ||
|
||
# References | ||
|
||
Roberts, Gareth O., and Jeffrey S. Rosenthal. "Examples of adaptive MCMC." | ||
Journal of Computational and Graphical Statistics 18.2 (2009): 349-367. | ||
""" | ||
mutable struct AdaptiveProposal{P} <: Proposal{P} | ||
proposal::P | ||
adaptor::Adaptor | ||
end | ||
|
||
function AdaptiveProposal(p; kwargs...) | ||
return AdaptiveProposal(p, Adaptor(; kwargs...)) | ||
end | ||
|
||
accepted!(p::AdaptiveProposal) = p.adaptor.accepted += 1 | ||
accepted!(p::Vector{<:AdaptiveProposal}) = map(accepted!, p) | ||
accepted!(p::NamedTuple{names}) where names = map(x->accepted!(getfield(p, x)), names) | ||
|
||
# this is defined because the first draw has no transition yet (I think) | ||
function propose(rng::Random.AbstractRNG, p::AdaptiveProposal, m::DensityModel) | ||
return rand(rng, p.proposal) | ||
end | ||
|
||
# the actual proposal happens here | ||
function propose( | ||
rng::Random.AbstractRNG, | ||
proposal::AdaptiveProposal{<:Union{Distribution,Proposal}}, | ||
model::DensityModel, | ||
t | ||
) | ||
consider_adaptation!(proposal) | ||
return t + rand(rng, proposal.proposal) | ||
end | ||
|
||
function q(proposal::AdaptiveProposal, t, t_cond) | ||
return logpdf(proposal, t - t_cond) | ||
end | ||
|
||
function consider_adaptation!(p) | ||
(p.adaptor.total % p.adaptor.tune == 0) && adapt!(p) | ||
p.adaptor.total += 1 | ||
end | ||
|
||
function adapt!(p::AdaptiveProposal) | ||
a = p.adaptor | ||
a.total == 0 && return | ||
δ = min(a.δmax, sqrt(a.tune / a.total)) # diminishing adaptation | ||
α = a.accepted / a.tune # acceptance ratio | ||
p_ = adapted(p.proposal, α > a.target ? δ : -δ, a.bound) | ||
a.accepted = 0 | ||
p.proposal = p_ | ||
end | ||
|
||
function adapted(d::Normal, δ, bound=Inf) | ||
_lσ = log(d.σ) + δ | ||
lσ = sign(_lσ) * min(bound, abs(_lσ)) | ||
return Normal(d.μ, exp(lσ)) | ||
end | ||
|
||
function adapted(d::Uniform, δ, bound=Inf) | ||
lσ = log(d.b) + δ | ||
σ = exp(sign(lσ) * min(bound, abs(lσ))) | ||
return Uniform(-σ, σ) | ||
end |
Original file line number | Diff line number | Diff line change |
---|---|---|
|
@@ -103,4 +103,4 @@ function q( | |
t_cond | ||
) | ||
return q(proposal(t_cond), t, t_cond) | ||
end | ||
end |
Original file line number | Diff line number | Diff line change |
---|---|---|
|
@@ -52,6 +52,32 @@ using Test | |
@test mean(chain2.μ) ≈ 0.0 atol=0.1 | ||
@test mean(chain2.σ) ≈ 1.0 atol=0.1 | ||
end | ||
|
||
@testset "Adaptive random walk" begin | ||
# Set up our sampler with initial parameters. | ||
spl1 = MetropolisHastings([AdaptiveProposal(Normal(0,.4)), AdaptiveProposal(Normal(0,1.2))]) | ||
spl2 = MetropolisHastings((μ=AdaptiveProposal(Normal(0,1.4)), σ=AdaptiveProposal(Normal(0,0.2)))) | ||
|
||
# Sample from the posterior. | ||
chain1 = sample(model, spl1, 100000; chain_type=StructArray, param_names=["μ", "σ"]) | ||
chain2 = sample(model, spl2, 100000; chain_type=StructArray, param_names=["μ", "σ"]) | ||
|
||
# chn_mean ≈ dist_mean atol=atol_v | ||
@test mean(chain1.μ) ≈ 0.0 atol=0.1 | ||
@test mean(chain1.σ) ≈ 1.0 atol=0.1 | ||
@test mean(chain2.μ) ≈ 0.0 atol=0.1 | ||
@test mean(chain2.σ) ≈ 1.0 atol=0.1 | ||
end | ||
|
||
@testset "Compare adaptive to simple random walk" begin | ||
data = rand(Normal(2., 1.), 500) | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. @arzwa This might be the problem - you redefine You could just rename the variable here but actually I think the better approach might be to "fix" the data in the model to avoid any such surprises in the future. I guess this can be achieved by defining density = let data = data
θ -> insupport(θ) ? sum(logpdf.(dist(θ), data)) : -Inf
end There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. But I haven't tested it, so make sure it actually fixes the problem 😄 There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Yes, I just saw this too, thanks. I'll check and push an updated test suite. (Actually, we could just as well test against the same data defined above in the test suite, but I find testing against a mean different from 0 a bit more reassuring since the sampler actually has to 'move' to somewhere from where it starts). |
||
m1 = DensityModel(x -> loglikelihood(Normal(x,1), data)) | ||
p1 = RandomWalkProposal(Normal()) | ||
p2 = AdaptiveProposal(Normal()) | ||
c1 = sample(m1, MetropolisHastings(p1), 10000; chain_type=Chains) | ||
c2 = sample(m1, MetropolisHastings(p2), 10000; chain_type=Chains) | ||
@test ess(c2).nt.ess > ess(c1).nt.ess | ||
end | ||
|
||
@testset "parallel sampling" begin | ||
spl1 = StaticMH([Normal(0,1), Normal(0, 1)]) | ||
|
@@ -105,4 +131,3 @@ using Test | |
|
||
@testset "EMCEE" begin include("emcee.jl") end | ||
end | ||
|
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Is this actually needed? It seems like the default for
Proposal
should cover this.There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Currently there seems to be no default
q
forProposal
. However it's also unnecessary because these adaptive proposals are always symmetric.