From e9fd60211c2f3721d379813977a2a30d863832f3 Mon Sep 17 00:00:00 2001 From: arzwa Date: Sun, 22 Nov 2020 16:59:31 +0100 Subject: [PATCH 01/31] adaptive proposal --- src/adaptive.jl | 88 ++++++++++++++++++++++++++++++++++++++++++++++++ test/runtests.jl | 18 +++++++++- 2 files changed, 105 insertions(+), 1 deletion(-) create mode 100644 src/adaptive.jl diff --git a/src/adaptive.jl b/src/adaptive.jl new file mode 100644 index 0000000..6773fca --- /dev/null +++ b/src/adaptive.jl @@ -0,0 +1,88 @@ +mutable struct Adaptor + accepted::Integer + total ::Integer + tune ::Integer # tuning interval + target ::Float64 # target acceptance rate + bound ::Float64 # bound on logσ of Gaussian kernel + δmax ::Float64 # maximum adaptation step +end + +Adaptor(; tune=25, target=0.44, bound=10., δmax=0.2) = + Adaptor(0, 0, tune, target, bound, δmax) + +""" + 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 +``` +""" +mutable struct AdaptiveProposal{P} <: Proposal{P} + proposal::P + adaptor ::Adaptor +end + +function AdaptiveProposal(p; kwargs...) + 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) +propose(rng::Random.AbstractRNG, p::AdaptiveProposal, m::DensityModel) = + rand(rng, p.proposal) + +# the actual proposal happens here +function propose( + rng::Random.AbstractRNG, + proposal::AdaptiveProposal{<:Union{Distribution,Proposal}}, + model::DensityModel, + t +) + consider_adaptation!(proposal) + t + rand(rng, proposal.proposal) +end + +function q(proposal::AdaptiveProposal, t, t_cond) + 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) + a = p.adaptor + a.total == 0 && return + δ = min(a.δmax, 1. /√(a.total/a.tune)) # 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σ = abs(lσ) > bound ? sign(lσ) * bound : lσ + Normal(d.μ, exp(lσ)) +end + +function adapted(d::Uniform, δ, bound=Inf) + lσ = log(d.b) + δ + σ = abs(lσ) > bound ? exp(sign(lσ) * bound) : exp(lσ) + Uniform(-σ, σ) +end + diff --git a/test/runtests.jl b/test/runtests.jl index 2e81db0..bbb6e31 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -11,7 +11,7 @@ using Test Random.seed!(1234) # Generate a set of data from the posterior we want to estimate. - data = rand(Normal(0, 1), 300) + data = rand(Normal(0., 1), 300) # Define the components of a basic model. insupport(θ) = θ[2] >= 0 @@ -52,6 +52,22 @@ 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 "parallel sampling" begin spl1 = StaticMH([Normal(0,1), Normal(0, 1)]) From 5f0ddfa406096909ec4ac94e3a1988e6ed8734d4 Mon Sep 17 00:00:00 2001 From: arzwa Date: Sun, 22 Nov 2020 17:50:54 +0100 Subject: [PATCH 02/31] Adaptive univariate proposals tests --- src/.adaptive.jl.swp | Bin 0 -> 12288 bytes src/AdvancedMH.jl | 1 + src/adaptive.jl | 7 ++++++- src/mh-core.jl | 4 ++++ src/proposal.jl | 2 +- test/runtests.jl | 9 +++++++++ 6 files changed, 21 insertions(+), 2 deletions(-) create mode 100644 src/.adaptive.jl.swp diff --git a/src/.adaptive.jl.swp b/src/.adaptive.jl.swp new file mode 100644 index 0000000000000000000000000000000000000000..75bfc9feb3bd3683b4ac88a8b83210590f548f00 GIT binary patch literal 12288 zcmeI2UuYaf9LJ~i|DUR~FFq*znqauq&E}FOZF4o1QZ+_REotqOlAF7oyW3`Oce^{6 zq?TxH@z0|u)W44^ErKEhK~WHVPnG&mEcj4-&=-r<2N5iOXJ+^Aa!FNu6ZXjG?sj&5 zf95;Cnc32-6c0^J(@v+z@VbMs3#UGt-MnF*z1q(hkDpoKs{J2c0k88!jA$2s*Y#QL zaL5(&LwnqL9#n-p{fKkeA4=lt5O;Yby?KG_zw?pb%N58K_zw!SWNmQgCbn&4Sg+Im zzOA%n^OOI*DR)6We(pXp`( zkt>iZkSmZYkSmZYkSmZYkSmZYkSmZYkSmZY@IO?5yNrE(4P)i&kUakXKl}IpPq#C6 z4txP#0~Ig;e!Y#cOWxj*mvMGcmtdO&w?dz5Db7muoe7xGh-isgWv$T8GLyYWAB6azzGn617HjI z<3`4Q1DC-HI0N1T?}8&>KX@EG2JQvF+`!m*@GUq8J_PT8<6t|u7JPL*V;_OH!Rz1@ zm;*E5ZmU_ z^)F+tBTrQx5>=^uF!e@wJsEhlFm4tG&z`jD>`6;Eo;^u4LA;U_d^FK+s)Vh;r3XkX zMuo(y2d4dA({52o-%x$EIljQ$59CjntDSR)SF0kD!bMw= zbXzp5bW~-F5qFH>_+zNb9m(UmkPvod8EFnG9k`>#wu}SxQyx05hm=(-5JdfC;;A33 zQOvPcbVzuDMbirk+*w&}@^(-q2&IM@nBiXhwGuMmo&VfZW2|Ff^vn$@*b~OwV&lT~kzNd$wJG4t31sS=uLN;?YNg50w(w~ z7St(Ii9(8b;1=SbZd2yJT_UDn*H164BFwhwCAARE&=r2E^q|0qj7e#vTclN;60`Y&vGgtjd{vlXgc%6G0SMKgotbiOClC z+hd{zk0NA64G{}GQQ_4!@U9VOZmt%Py0-HyNt>Fh$84dqRLx| z(D~K8IFaPWd;zH?6oiCRK&hrqj-yywg(>v&f;rS*wmfxx?&D429&AOvI2Pb*vSiVM z8Wia(wSst5SRCBlF~UNvdBVp!#(`s5w&`iDcB}T+Fveu;NzbZmaPV(@wjH8UdW?H* zlu2hvWQJzFw5z(;NTsLS!nm-D<`%HSB#z@)T_wKOq`PBRVRJQzFs8KJu+r6xo4MsE z(EY60lDy&zLd~_R(gen_myP<^{w2DAll?oj!&-n-Mh>0Cv*^5|qnjR&jKZjRtunjdUESK4y^tn*appdG+?Puf{HG!;-7<5(8PK-r<3~@lwAJi4+0%tE t2KvcE1)B<_s6$ni9wbp^=b 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 @@ -64,7 +69,7 @@ function consider_adaptation!(p) p.adaptor.total += 1 end -function adapt!(p) +function adapt!(p::AdaptiveProposal) a = p.adaptor a.total == 0 && return δ = min(a.δmax, 1. /√(a.total/a.tune)) # diminishing adaptation diff --git a/src/mh-core.jl b/src/mh-core.jl index 33832a2..39560b1 100644 --- a/src/mh-core.jl +++ b/src/mh-core.jl @@ -208,8 +208,12 @@ function AbstractMCMC.step( # Decide whether to return the previous params or the new one. if -Random.randexp(rng) < logα + accepted!(spl.proposal) return params, params else return params_prev, params_prev end end + +function accepted!(p::P) where P<:Proposal end + diff --git a/src/proposal.jl b/src/proposal.jl index 286ad94..86d9462 100644 --- a/src/proposal.jl +++ b/src/proposal.jl @@ -103,4 +103,4 @@ function q( t_cond ) return q(proposal(t_cond), t, t_cond) -end \ No newline at end of file +end diff --git a/test/runtests.jl b/test/runtests.jl index bbb6e31..21a2607 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -69,6 +69,15 @@ using Test @test mean(chain2.σ) ≈ 1.0 atol=0.1 end + @testset "Compare adaptive to simple random walk" begin + data = rand(Normal(2., 1.), 500) + 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) + end + @testset "parallel sampling" begin spl1 = StaticMH([Normal(0,1), Normal(0, 1)]) From f42b78475642c176d370e5809b16ff2e2a5b6631 Mon Sep 17 00:00:00 2001 From: arzwa Date: Sun, 22 Nov 2020 17:54:22 +0100 Subject: [PATCH 03/31] swp files from vim... --- src/.adaptive.jl.swp | Bin 12288 -> 0 bytes 1 file changed, 0 insertions(+), 0 deletions(-) delete mode 100644 src/.adaptive.jl.swp diff --git a/src/.adaptive.jl.swp b/src/.adaptive.jl.swp deleted file mode 100644 index 75bfc9feb3bd3683b4ac88a8b83210590f548f00..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 12288 zcmeI2UuYaf9LJ~i|DUR~FFq*znqauq&E}FOZF4o1QZ+_REotqOlAF7oyW3`Oce^{6 zq?TxH@z0|u)W44^ErKEhK~WHVPnG&mEcj4-&=-r<2N5iOXJ+^Aa!FNu6ZXjG?sj&5 zf95;Cnc32-6c0^J(@v+z@VbMs3#UGt-MnF*z1q(hkDpoKs{J2c0k88!jA$2s*Y#QL zaL5(&LwnqL9#n-p{fKkeA4=lt5O;Yby?KG_zw?pb%N58K_zw!SWNmQgCbn&4Sg+Im zzOA%n^OOI*DR)6We(pXp`( zkt>iZkSmZYkSmZYkSmZYkSmZYkSmZYkSmZY@IO?5yNrE(4P)i&kUakXKl}IpPq#C6 z4txP#0~Ig;e!Y#cOWxj*mvMGcmtdO&w?dz5Db7muoe7xGh-isgWv$T8GLyYWAB6azzGn617HjI z<3`4Q1DC-HI0N1T?}8&>KX@EG2JQvF+`!m*@GUq8J_PT8<6t|u7JPL*V;_OH!Rz1@ zm;*E5ZmU_ z^)F+tBTrQx5>=^uF!e@wJsEhlFm4tG&z`jD>`6;Eo;^u4LA;U_d^FK+s)Vh;r3XkX zMuo(y2d4dA({52o-%x$EIljQ$59CjntDSR)SF0kD!bMw= zbXzp5bW~-F5qFH>_+zNb9m(UmkPvod8EFnG9k`>#wu}SxQyx05hm=(-5JdfC;;A33 zQOvPcbVzuDMbirk+*w&}@^(-q2&IM@nBiXhwGuMmo&VfZW2|Ff^vn$@*b~OwV&lT~kzNd$wJG4t31sS=uLN;?YNg50w(w~ z7St(Ii9(8b;1=SbZd2yJT_UDn*H164BFwhwCAARE&=r2E^q|0qj7e#vTclN;60`Y&vGgtjd{vlXgc%6G0SMKgotbiOClC z+hd{zk0NA64G{}GQQ_4!@U9VOZmt%Py0-HyNt>Fh$84dqRLx| z(D~K8IFaPWd;zH?6oiCRK&hrqj-yywg(>v&f;rS*wmfxx?&D429&AOvI2Pb*vSiVM z8Wia(wSst5SRCBlF~UNvdBVp!#(`s5w&`iDcB}T+Fveu;NzbZmaPV(@wjH8UdW?H* zlu2hvWQJzFw5z(;NTsLS!nm-D<`%HSB#z@)T_wKOq`PBRVRJQzFs8KJu+r6xo4MsE z(EY60lDy&zLd~_R(gen_myP<^{w2DAll?oj!&-n-Mh>0Cv*^5|qnjR&jKZjRtunjdUESK4y^tn*appdG+?Puf{HG!;-7<5(8PK-r<3~@lwAJi4+0%tE t2KvcE1)B<_s6$ni9wbp^=b Date: Sun, 22 Nov 2020 17:58:40 +0100 Subject: [PATCH 04/31] test ESS --- test/runtests.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/test/runtests.jl b/test/runtests.jl index 21a2607..d129805 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -76,6 +76,7 @@ using Test 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 From d8989aa9fa3896ee059d1867c3646abcf6905db8 Mon Sep 17 00:00:00 2001 From: arzwa Date: Sun, 22 Nov 2020 19:08:33 +0100 Subject: [PATCH 05/31] export AdaptiveProposal --- src/AdvancedMH.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/AdvancedMH.jl b/src/AdvancedMH.jl index c13c8e6..f2a3918 100644 --- a/src/AdvancedMH.jl +++ b/src/AdvancedMH.jl @@ -9,7 +9,7 @@ import Random # Exports export MetropolisHastings, DensityModel, RWMH, StaticMH, StaticProposal, - RandomWalkProposal, Ensemble, StretchProposal + RandomWalkProposal, Ensemble, StretchProposal, AdaptiveProposal # Reexports export sample, MCMCThreads, MCMCDistributed From 565f12a9d7e04d4900a79e2ec368510b2256e67a Mon Sep 17 00:00:00 2001 From: Arthur Zwaenepoel Date: Sun, 22 Nov 2020 19:58:41 +0100 Subject: [PATCH 06/31] Update src/adaptive.jl Co-authored-by: Cameron Pfiffer --- src/adaptive.jl | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/src/adaptive.jl b/src/adaptive.jl index 3e37bd1..b5018ec 100644 --- a/src/adaptive.jl +++ b/src/adaptive.jl @@ -1,10 +1,10 @@ mutable struct Adaptor accepted::Integer - total ::Integer - tune ::Integer # tuning interval - target ::Float64 # target acceptance rate - bound ::Float64 # bound on logσ of Gaussian kernel - δmax ::Float64 # maximum adaptation step + total::Integer + tune::Integer # tuning interval + target::Float64 # target acceptance rate + bound::Float64 # bound on logσ of Gaussian kernel + δmax::Float64 # maximum adaptation step end Adaptor(; tune=25, target=0.44, bound=10., δmax=0.2) = @@ -90,4 +90,3 @@ function adapted(d::Uniform, δ, bound=Inf) σ = abs(lσ) > bound ? exp(sign(lσ) * bound) : exp(lσ) Uniform(-σ, σ) end - From cc16195968faded0487517beb101ba545c2d81bd Mon Sep 17 00:00:00 2001 From: Arthur Zwaenepoel Date: Sun, 22 Nov 2020 19:59:03 +0100 Subject: [PATCH 07/31] Update src/adaptive.jl Co-authored-by: Cameron Pfiffer --- src/adaptive.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/adaptive.jl b/src/adaptive.jl index b5018ec..8f916fd 100644 --- a/src/adaptive.jl +++ b/src/adaptive.jl @@ -34,7 +34,7 @@ Journal of Computational and Graphical Statistics 18.2 (2009): 349-367. """ mutable struct AdaptiveProposal{P} <: Proposal{P} proposal::P - adaptor ::Adaptor + adaptor::Adaptor end function AdaptiveProposal(p; kwargs...) From 802ec67b6a2875ac264377dc5d9977a52c64f01d Mon Sep 17 00:00:00 2001 From: Arthur Zwaenepoel Date: Sun, 22 Nov 2020 19:59:12 +0100 Subject: [PATCH 08/31] Update src/adaptive.jl Co-authored-by: Cameron Pfiffer --- src/adaptive.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/adaptive.jl b/src/adaptive.jl index 8f916fd..2b62524 100644 --- a/src/adaptive.jl +++ b/src/adaptive.jl @@ -46,8 +46,9 @@ 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) -propose(rng::Random.AbstractRNG, p::AdaptiveProposal, m::DensityModel) = - rand(rng, p.proposal) +function propose(rng::Random.AbstractRNG, p::AdaptiveProposal, m::DensityModel) + return rand(rng, p.proposal) +end # the actual proposal happens here function propose( From c7623c48d0cd94eff031c516a93e16e58ac355e9 Mon Sep 17 00:00:00 2001 From: Arthur Zwaenepoel Date: Sun, 22 Nov 2020 19:59:26 +0100 Subject: [PATCH 09/31] Update src/adaptive.jl Co-authored-by: Cameron Pfiffer --- src/adaptive.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/adaptive.jl b/src/adaptive.jl index 2b62524..010545f 100644 --- a/src/adaptive.jl +++ b/src/adaptive.jl @@ -7,8 +7,9 @@ mutable struct Adaptor δmax::Float64 # maximum adaptation step end -Adaptor(; tune=25, target=0.44, bound=10., δmax=0.2) = - Adaptor(0, 0, tune, target, bound, δmax) +function Adaptor(; tune=25, target=0.44, bound=10., δmax=0.2) + return Adaptor(0, 0, tune, target, bound, δmax) +end """ AdaptiveProposal{P} From a33937eb50b51d781dde16fddccd84ba5d4b4612 Mon Sep 17 00:00:00 2001 From: Arthur Zwaenepoel Date: Sun, 22 Nov 2020 19:59:34 +0100 Subject: [PATCH 10/31] Update src/adaptive.jl Co-authored-by: Cameron Pfiffer --- src/adaptive.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/adaptive.jl b/src/adaptive.jl index 010545f..7fa2cf2 100644 --- a/src/adaptive.jl +++ b/src/adaptive.jl @@ -39,7 +39,7 @@ mutable struct AdaptiveProposal{P} <: Proposal{P} end function AdaptiveProposal(p; kwargs...) - AdaptiveProposal(p, Adaptor(; kwargs...)) + return AdaptiveProposal(p, Adaptor(; kwargs...)) end accepted!(p::AdaptiveProposal) = p.adaptor.accepted += 1 From 4007bd0c9865832e8b6f526dc998fd332cce3f83 Mon Sep 17 00:00:00 2001 From: Arthur Zwaenepoel Date: Sun, 22 Nov 2020 19:59:44 +0100 Subject: [PATCH 11/31] Update src/adaptive.jl Co-authored-by: Cameron Pfiffer --- src/adaptive.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/adaptive.jl b/src/adaptive.jl index 7fa2cf2..0e9dbb0 100644 --- a/src/adaptive.jl +++ b/src/adaptive.jl @@ -59,7 +59,7 @@ function propose( t ) consider_adaptation!(proposal) - t + rand(rng, proposal.proposal) + return t + rand(rng, proposal.proposal) end function q(proposal::AdaptiveProposal, t, t_cond) From 71e010b89dd2ecc48f7b529e4eac61fdf6a43d49 Mon Sep 17 00:00:00 2001 From: Arthur Zwaenepoel Date: Sun, 22 Nov 2020 19:59:52 +0100 Subject: [PATCH 12/31] Update src/adaptive.jl Co-authored-by: Cameron Pfiffer --- src/adaptive.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/adaptive.jl b/src/adaptive.jl index 0e9dbb0..4fa2946 100644 --- a/src/adaptive.jl +++ b/src/adaptive.jl @@ -63,7 +63,7 @@ function propose( end function q(proposal::AdaptiveProposal, t, t_cond) - logpdf(proposal, t - t_cond) + return logpdf(proposal, t - t_cond) end function consider_adaptation!(p) From 2d59ede8f1f9b5d511a8f50ffb7aa443feb2f257 Mon Sep 17 00:00:00 2001 From: Arthur Zwaenepoel Date: Sun, 22 Nov 2020 20:00:03 +0100 Subject: [PATCH 13/31] Update src/adaptive.jl Co-authored-by: Cameron Pfiffer --- src/adaptive.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/adaptive.jl b/src/adaptive.jl index 4fa2946..c24f0da 100644 --- a/src/adaptive.jl +++ b/src/adaptive.jl @@ -84,7 +84,7 @@ end function adapted(d::Normal, δ, bound=Inf) lσ = log(d.σ) + δ lσ = abs(lσ) > bound ? sign(lσ) * bound : lσ - Normal(d.μ, exp(lσ)) + return Normal(d.μ, exp(lσ)) end function adapted(d::Uniform, δ, bound=Inf) From 279aea747f3c61938e98631aadd4f8c699bc5b73 Mon Sep 17 00:00:00 2001 From: Arthur Zwaenepoel Date: Sun, 22 Nov 2020 20:00:13 +0100 Subject: [PATCH 14/31] Update src/adaptive.jl Co-authored-by: Cameron Pfiffer --- src/adaptive.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/adaptive.jl b/src/adaptive.jl index c24f0da..440420a 100644 --- a/src/adaptive.jl +++ b/src/adaptive.jl @@ -90,5 +90,5 @@ end function adapted(d::Uniform, δ, bound=Inf) lσ = log(d.b) + δ σ = abs(lσ) > bound ? exp(sign(lσ) * bound) : exp(lσ) - Uniform(-σ, σ) + return Uniform(-σ, σ) end From 93f17c579f9280c5df7f255f28bcce240e9cbae2 Mon Sep 17 00:00:00 2001 From: Arthur Zwaenepoel Date: Mon, 23 Nov 2020 08:12:29 +0100 Subject: [PATCH 15/31] Update src/adaptive.jl Co-authored-by: David Widmann --- src/adaptive.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/adaptive.jl b/src/adaptive.jl index 440420a..cf41e0e 100644 --- a/src/adaptive.jl +++ b/src/adaptive.jl @@ -89,6 +89,6 @@ end function adapted(d::Uniform, δ, bound=Inf) lσ = log(d.b) + δ - σ = abs(lσ) > bound ? exp(sign(lσ) * bound) : exp(lσ) + σ = exp(sign(lσ) * min(bound, abs(lσ))) return Uniform(-σ, σ) end From 16715e1bd6b6e3818e92bc9e257861697fb94236 Mon Sep 17 00:00:00 2001 From: Arthur Zwaenepoel Date: Mon, 23 Nov 2020 08:13:22 +0100 Subject: [PATCH 16/31] Update src/adaptive.jl Co-authored-by: David Widmann --- src/adaptive.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/adaptive.jl b/src/adaptive.jl index cf41e0e..330ad4a 100644 --- a/src/adaptive.jl +++ b/src/adaptive.jl @@ -82,8 +82,8 @@ function adapt!(p::AdaptiveProposal) end function adapted(d::Normal, δ, bound=Inf) - lσ = log(d.σ) + δ - lσ = abs(lσ) > bound ? sign(lσ) * bound : lσ + _lσ = log(d.σ) + δ + lσ = sign(_lσ) * min(bound, abs(_lσ)) return Normal(d.μ, exp(lσ)) end From 046c21bfcaa1fe2b56bafa393aaa8b092e0aa8cf Mon Sep 17 00:00:00 2001 From: Arthur Zwaenepoel Date: Mon, 23 Nov 2020 08:13:34 +0100 Subject: [PATCH 17/31] Update src/adaptive.jl Co-authored-by: David Widmann --- src/adaptive.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/adaptive.jl b/src/adaptive.jl index 330ad4a..66a6757 100644 --- a/src/adaptive.jl +++ b/src/adaptive.jl @@ -74,7 +74,7 @@ end function adapt!(p::AdaptiveProposal) a = p.adaptor a.total == 0 && return - δ = min(a.δmax, 1. /√(a.total/a.tune)) # diminishing adaptation + δ = 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 From b91fcc03ddea53d7aebedc2770e38720c9963757 Mon Sep 17 00:00:00 2001 From: arzwa Date: Mon, 23 Nov 2020 08:35:33 +0100 Subject: [PATCH 18/31] test update --- test/runtests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index d129805..b47fef4 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -75,7 +75,7 @@ using Test p1 = RandomWalkProposal(Normal()) p2 = AdaptiveProposal(Normal()) c1 = sample(m1, MetropolisHastings(p1), 10000; chain_type=Chains) - c2 = sample(m1, MetropolisHastings(p2), 10000; chain_type=Chains) + c2 = sample(m1, MetropolisHastings(p2), 10000; chain_type=Chains) @test ess(c2).nt.ess > ess(c1).nt.ess end From 8fb104829bb657a443690459e2fa4d0cdab16a27 Mon Sep 17 00:00:00 2001 From: arzwa Date: Mon, 23 Nov 2020 08:51:19 +0100 Subject: [PATCH 19/31] Adaptor docstring --- src/adaptive.jl | 19 ++++++++++++++++--- 1 file changed, 16 insertions(+), 3 deletions(-) diff --git a/src/adaptive.jl b/src/adaptive.jl index 66a6757..795c8ed 100644 --- a/src/adaptive.jl +++ b/src/adaptive.jl @@ -1,7 +1,20 @@ +""" + 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::Integer - total::Integer - tune::Integer # tuning interval + 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 From 407167599f64f221e478fa510ad218346497de04 Mon Sep 17 00:00:00 2001 From: Arthur Zwaenepoel Date: Mon, 23 Nov 2020 11:45:18 +0100 Subject: [PATCH 20/31] Update test/runtests.jl Co-authored-by: David Widmann --- test/runtests.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index b47fef4..d0fbd46 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -11,7 +11,7 @@ using Test Random.seed!(1234) # Generate a set of data from the posterior we want to estimate. - data = rand(Normal(0., 1), 300) + data = rand(Normal(0, 1), 300) # Define the components of a basic model. insupport(θ) = θ[2] >= 0 @@ -131,4 +131,3 @@ using Test @testset "EMCEE" begin include("emcee.jl") end end - From a63262d651c3e0a25f49d2c06249f699eadc5bf2 Mon Sep 17 00:00:00 2001 From: arzwa Date: Mon, 23 Nov 2020 12:10:58 +0100 Subject: [PATCH 21/31] updated tests --- test/emcee.jl | 4 ++-- test/runtests.jl | 50 +++++++++++++++++++++++++++++------------------- 2 files changed, 32 insertions(+), 22 deletions(-) diff --git a/test/emcee.jl b/test/emcee.jl index 0803ed3..039e9af 100644 --- a/test/emcee.jl +++ b/test/emcee.jl @@ -18,7 +18,7 @@ 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) + param_names = ["s", "m"], chain_type = Chains, progress = false) @test mean(chain["s"]) ≈ 49/24 atol=0.1 @test mean(chain["m"]) ≈ 7/6 atol=0.1 @@ -43,7 +43,7 @@ Random.seed!(100) sampler = Ensemble(1_000, StretchProposal(MvNormal(2, 1))) chain = sample(model, sampler, 1_000; - param_names = ["logs", "m"], chain_type = Chains) + param_names = ["logs", "m"], chain_type = Chains, progress = false) @test mean(exp, chain["logs"]) ≈ 49/24 atol=0.1 @test mean(chain["m"]) ≈ 7/6 atol=0.1 diff --git a/test/runtests.jl b/test/runtests.jl index b47fef4..bd2d30d 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -16,7 +16,11 @@ using Test # Define the components of a basic model. insupport(θ) = θ[2] >= 0 dist(θ) = Normal(θ[1], θ[2]) - density(θ) = insupport(θ) ? sum(logpdf.(dist(θ), data)) : -Inf + + # using the let statement prevents surprises when data gets redefined + density = let data = data + θ -> insupport(θ) ? sum(logpdf.(dist(θ), data)) : -Inf + end # Construct a DensityModel. model = DensityModel(density) @@ -27,8 +31,9 @@ using Test spl2 = StaticMH(MvNormal([0.0, 0.0], 1)) # Sample from the posterior. - chain1 = sample(model, spl1, 100000; chain_type=StructArray, param_names=["μ", "σ"]) - chain2 = sample(model, spl2, 100000; chain_type=StructArray, param_names=["μ", "σ"]) + kwargs = (progress=false, chain_type=StructArray, param_names=["μ", "σ"]) + chain1 = sample(model, spl1, 100000; kwargs...) + chain2 = sample(model, spl2, 100000; kwargs...) # chn_mean ≈ dist_mean atol=atol_v @test mean(chain1.μ) ≈ 0.0 atol=0.1 @@ -43,8 +48,9 @@ using Test spl2 = RWMH(MvNormal([0.0, 0.0], 1)) # Sample from the posterior. - chain1 = sample(model, spl1, 100000; chain_type=StructArray, param_names=["μ", "σ"]) - chain2 = sample(model, spl2, 100000; chain_type=StructArray, param_names=["μ", "σ"]) + kwargs = (progress=false, chain_type=StructArray, param_names=["μ", "σ"]) + chain1 = sample(model, spl1, 100000; kwargs...) + chain2 = sample(model, spl2, 100000; kwargs...) # chn_mean ≈ dist_mean atol=atol_v @test mean(chain1.μ) ≈ 0.0 atol=0.1 @@ -55,12 +61,15 @@ using Test @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)))) + p1 = [AdaptiveProposal(Normal(0,.4)), AdaptiveProposal(Normal(0,1.2))] + p2 = (μ=AdaptiveProposal(Normal(0,1.4)), σ=AdaptiveProposal(Normal(0,0.2))) + spl1 = MetropolisHastings(p1) + spl2 = MetropolisHastings(p2) # Sample from the posterior. - chain1 = sample(model, spl1, 100000; chain_type=StructArray, param_names=["μ", "σ"]) - chain2 = sample(model, spl2, 100000; chain_type=StructArray, param_names=["μ", "σ"]) + kwargs = (progress=false, chain_type=StructArray, param_names=["μ", "σ"]) + chain1 = sample(model, spl1, 100000; kwargs...) + chain2 = sample(model, spl2, 100000; kwargs...) # chn_mean ≈ dist_mean atol=atol_v @test mean(chain1.μ) ≈ 0.0 atol=0.1 @@ -74,22 +83,22 @@ using Test 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) + kwargs = (progress=false, chain_type=Chains) + c1 = sample(m1, MetropolisHastings(p1), 10000; kwargs...) + c2 = sample(m1, MetropolisHastings(p2), 10000; kwargs...) @test ess(c2).nt.ess > ess(c1).nt.ess end @testset "parallel sampling" begin spl1 = StaticMH([Normal(0,1), Normal(0, 1)]) - chain1 = sample(model, spl1, MCMCDistributed(), 10000, 4; - param_names=["μ", "σ"], chain_type=Chains) + kwargs = (progress=false, chain_type=Chains, param_names=["μ", "σ"]) + chain1 = sample(model, spl1, MCMCDistributed(), 10000, 4; kwargs...) @test mean(chain1["μ"]) ≈ 0.0 atol=0.1 @test mean(chain1["σ"]) ≈ 1.0 atol=0.1 if VERSION >= v"1.3" - chain2 = sample(model, spl1, MCMCThreads(), 10000, 4; - param_names=["μ", "σ"], chain_type=Chains) + chain2 = sample(model, spl1, MCMCThreads(), 10000, 4; kwargs...) @test mean(chain2["μ"]) ≈ 0.0 atol=0.1 @test mean(chain2["σ"]) ≈ 1.0 atol=0.1 end @@ -106,10 +115,11 @@ using Test p3 = (a=StaticProposal(Normal(0,1)), b=StaticProposal(InverseGamma(2,3))) p4 = StaticProposal((x=1.0) -> Normal(x, 1)) - c1 = sample(m1, MetropolisHastings(p1), 100; chain_type=Vector{NamedTuple}) - c2 = sample(m2, MetropolisHastings(p2), 100; chain_type=Vector{NamedTuple}) - c3 = sample(m3, MetropolisHastings(p3), 100; chain_type=Vector{NamedTuple}) - c4 = sample(m4, MetropolisHastings(p4), 100; chain_type=Vector{NamedTuple}) + kwargs = (chain_type=Vector{NamedTuple}, progress=false) + c1 = sample(m1, MetropolisHastings(p1), 100; kwargs...) + c2 = sample(m2, MetropolisHastings(p2), 100; kwargs...) + c3 = sample(m3, MetropolisHastings(p3), 100; kwargs...) + c4 = sample(m4, MetropolisHastings(p4), 100; kwargs...) @test keys(c1[1]) == (:param_1, :lp) @test keys(c2[1]) == (:param_1, :param_2, :lp) @@ -124,7 +134,7 @@ using Test val = [0.4, 1.2] # Sample from the posterior. - chain1 = sample(model, spl1, 10, init_params = val) + chain1 = sample(model, spl1, 10, init_params = val, progress=false) @test chain1[1].params == val end From afd3ed196fa16822c34d77ec069102aaa212544a Mon Sep 17 00:00:00 2001 From: arzwa Date: Mon, 23 Nov 2020 12:12:21 +0100 Subject: [PATCH 22/31] updated tests --- test/runtests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index bd2d30d..560ebfe 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -17,7 +17,7 @@ using Test insupport(θ) = θ[2] >= 0 dist(θ) = Normal(θ[1], θ[2]) - # using the let statement prevents surprises when data gets redefined + # using `let` prevents surprises when data is redefined in some testset density = let data = data θ -> insupport(θ) ? sum(logpdf.(dist(θ), data)) : -Inf end From f66f647e6369aa67f90df5274d96494a6fb1a16d Mon Sep 17 00:00:00 2001 From: arzwa Date: Thu, 10 Dec 2020 21:12:51 +0100 Subject: [PATCH 23/31] adaptive MvNormal --- Project.toml | 2 + src/AdvancedMH.jl | 6 ++- src/adaptivemvnormal.jl | 82 +++++++++++++++++++++++++++++++++++++++++ test/Project.toml | 4 ++ test/runtests.jl | 22 +++++++++++ 5 files changed, 115 insertions(+), 1 deletion(-) create mode 100644 src/adaptivemvnormal.jl create mode 100644 test/Project.toml diff --git a/Project.toml b/Project.toml index 6ba3217..5943307 100644 --- a/Project.toml +++ b/Project.toml @@ -5,6 +5,8 @@ version = "0.5.4" [deps] AbstractMCMC = "80f14c24-f653-4e6a-9b94-39d6b0f70001" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +PDMats = "90014a1f-27ba-587c-ab20-58faa44d9150" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Requires = "ae029012-a4dd-5104-9daa-d747884805df" diff --git a/src/AdvancedMH.jl b/src/AdvancedMH.jl index f2a3918..8014be5 100644 --- a/src/AdvancedMH.jl +++ b/src/AdvancedMH.jl @@ -4,12 +4,15 @@ module AdvancedMH using AbstractMCMC using Distributions using Requires +using LinearAlgebra +using PDMats import Random # Exports export MetropolisHastings, DensityModel, RWMH, StaticMH, StaticProposal, - RandomWalkProposal, Ensemble, StretchProposal, AdaptiveProposal + RandomWalkProposal, Ensemble, StretchProposal, AdaptiveProposal, + AdaptiveMvNormal # Reexports export sample, MCMCThreads, MCMCDistributed @@ -105,5 +108,6 @@ include("proposal.jl") include("mh-core.jl") include("emcee.jl") include("adaptive.jl") +include("adaptivemvnormal.jl") end # module AdvancedMH diff --git a/src/adaptivemvnormal.jl b/src/adaptivemvnormal.jl new file mode 100644 index 0000000..d02648b --- /dev/null +++ b/src/adaptivemvnormal.jl @@ -0,0 +1,82 @@ +""" + AdaptiveMvNormal(constant_component; σ=2.38, β=0.05) + +An adaptive multivariate normal proposal. More precisely this uses a +two-component mixture (with weights `β` and `1 - β`) of multivariate normal +proposal distributions, where the low-weight component is constant (the +`constant` field and the high-weight component has its covariance matrix +adapted to the covariance structure of the target (the `adaptive` field). +""" +mutable struct AdaptiveMvNormal{T1,T2,V} <: Proposal{T1} + d::Int # dimensionality + n::Int # iteration + β::Float64 # constant component mixture weight + σ::Float64 # scale factor for adapted covariance matrix + constant::T1 + adaptive::T2 + Ex::Vector{V} # rolling mean vector + EX::Matrix{V} # scatter matrix of previous draws +end + +function AdaptiveMvNormal(dist; σ=2.38, β=0.05) + n = length(dist) + # the `adaptive` MvNormal part must have PDMat covariance matrix, not + # ScalMat or PDDiagMat (as the `constant` part will have) + adaptive = MvNormal(PDMat(Symmetric(dist.Σ))) + AdaptiveMvNormal(n, -1, β, σ, dist, adaptive, zeros(n), zeros(n,n)) +end + +""" + adapt!(p::AdaptiveMvNormal, x::AbstractVector) + +Adaptation for the adaptive multivariate normal mixture proposal as described +in Haario et al. and Roberts & Rosenthal (2009). The code for this routine is +largely taken from `Mamba.jl`. + +# References + +Roberts, Gareth O., and Jeffrey S. Rosenthal. "Examples of adaptive MCMC." +Journal of Computational and Graphical Statistics 18.2 (2009): 349-367. +""" +function adapt!(p::AdaptiveMvNormal, x::AbstractVector) + p.n += 1 + # adapt mean vector and scatter matrix + f = p.n / (p.n + 1.0) + p.Ex = f * p.Ex + (1.0 - f) * x + p.EX = f * p.EX + (1.0 - f) * x * x' + # compute adapted covariance matrix + Σ = (p.σ^2 / (p.d * f)) * (p.EX - p.Ex * p.Ex') + #F = cholesky(Hermitian(Σ), Val{true}(), check=false) + F = cholesky(Hermitian(Σ), check=false) + if rank(F.L) == p.d + #p.adaptive = MvNormal(PDMat(Symmetric(F.P * F.L))) + p.adaptive = MvNormal(PDMat(Σ, F)) + end +end + +function propose(rng::Random.AbstractRNG, p::AdaptiveMvNormal) + return if p.n > 2p.d + p.β * rand(rng, p.constant) + (1.0 - p.β) * rand(rng, p.adaptive) + else + rand(rng, p.constant) + end +end + +function propose(rng::Random.AbstractRNG, p::AdaptiveMvNormal, m::DensityModel) + return propose(rng, p) +end + +function propose( + rng::Random.AbstractRNG, + proposal::AdaptiveMvNormal, + model::DensityModel, + t +) + adapt!(proposal, t) + return t + propose(rng, proposal) +end + +# this is always symmetric, so ignore +function q(proposal::AdaptiveMvNormal, t, t_cond) + return 0. +end diff --git a/test/Project.toml b/test/Project.toml new file mode 100644 index 0000000..5b194a1 --- /dev/null +++ b/test/Project.toml @@ -0,0 +1,4 @@ +[deps] +AdvancedMH = "5b7e9947-ddc0-4b3f-9b55-0d8042f74170" +MCMCChains = "c7f686f2-ff18-58e9-bc7b-31028e88f75d" +StructArrays = "09ab397b-f2b6-538f-b94a-2f83cf4a842a" diff --git a/test/runtests.jl b/test/runtests.jl index d16a7ba..5d2c90d 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -88,6 +88,28 @@ using Test c2 = sample(m1, MetropolisHastings(p2), 10000; kwargs...) @test ess(c2).nt.ess > ess(c1).nt.ess end + + @testset "Adaptive MvNormal mixture" begin + p1 = AdaptiveMvNormal(MvNormal(2, 0.1)) + spl1 = MetropolisHastings(p1) + kwargs = (progress=false, chain_type=StructArray, param_names=["μ", "σ"]) + chain1 = sample(model, spl1, 100000; kwargs...) + @test mean(chain1.μ) ≈ 0.0 atol=0.1 + @test mean(chain1.σ) ≈ 1.0 atol=0.1 + end + + @testset "Adaptive MvNormal mixture ESS" begin + Random.seed!(12) + d = 25 + M = randn(d,d) + Σ = M*M' + m = DensityModel(x -> logpdf(MvNormal(Σ), x)) + p = AdaptiveMvNormal(MvNormal(d, 1.)) + kwargs = (progress=false, chain_type=Chains) + c1 = sample(m, MetropolisHastings(p), 10000; kwargs...) + c2 = sample(m, RWMH(MvNormal(zeros(d), 1)), 10000; kwargs...) + @test sum(ess(c1).nt.ess .> ess(c2).nt.ess) > 20 + end @testset "parallel sampling" begin spl1 = StaticMH([Normal(0,1), Normal(0, 1)]) From e5ad0410580cccd475467fd54a0b336f0cca5629 Mon Sep 17 00:00:00 2001 From: arzwa Date: Fri, 11 Dec 2020 17:48:33 +0100 Subject: [PATCH 24/31] AdaptiveMvNormal with tests --- src/adaptive.jl | 3 +++ src/adaptivemvnormal.jl | 54 +++++++++++++++++++++-------------------- test/runtests.jl | 15 ++++++------ 3 files changed, 39 insertions(+), 33 deletions(-) diff --git a/src/adaptive.jl b/src/adaptive.jl index 795c8ed..536e208 100644 --- a/src/adaptive.jl +++ b/src/adaptive.jl @@ -55,6 +55,9 @@ function AdaptiveProposal(p; kwargs...) return AdaptiveProposal(p, Adaptor(; kwargs...)) end +# Adaptive proposals are only defined for symmetric proposal distributions +is_symmetric_proposal(::AdaptiveProposal) = true + 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) diff --git a/src/adaptivemvnormal.jl b/src/adaptivemvnormal.jl index d02648b..bf0372d 100644 --- a/src/adaptivemvnormal.jl +++ b/src/adaptivemvnormal.jl @@ -1,11 +1,22 @@ """ - AdaptiveMvNormal(constant_component; σ=2.38, β=0.05) + AdaptiveMvNormal(constant_component::MvNormal; σ=2.38, β=0.05) -An adaptive multivariate normal proposal. More precisely this uses a -two-component mixture (with weights `β` and `1 - β`) of multivariate normal -proposal distributions, where the low-weight component is constant (the -`constant` field and the high-weight component has its covariance matrix -adapted to the covariance structure of the target (the `adaptive` field). +Adaptive multivariate normal mixture proposal as described in Haario et al. and +Roberts & Rosenthal (2009). Uses a two-component mixture of MvNormal +distributions. One of the components (with mixture weight `β`) remains +constant, while the other component is adapted to the target covariance +structure. The proposal is initialized by providing the constant component to +the constructor. + +`σ` is the scale factor for the covariance matrix, where 2.38 is supposedly +optimal in a high-dimensional context according to Roberts & Rosenthal. + +# References + +- Haario, Heikki, Eero Saksman, and Johanna Tamminen. + "An adaptive Metropolis algorithm." Bernoulli 7.2 (2001): 223-242. +- Roberts, Gareth O., and Jeffrey S. Rosenthal. "Examples of adaptive MCMC." + Journal of Computational and Graphical Statistics 18.2 (2009): 349-367. """ mutable struct AdaptiveMvNormal{T1,T2,V} <: Proposal{T1} d::Int # dimensionality @@ -18,25 +29,21 @@ mutable struct AdaptiveMvNormal{T1,T2,V} <: Proposal{T1} EX::Matrix{V} # scatter matrix of previous draws end -function AdaptiveMvNormal(dist; σ=2.38, β=0.05) +function AdaptiveMvNormal(dist::MvNormal; σ=2.38, β=0.05) n = length(dist) - # the `adaptive` MvNormal part must have PDMat covariance matrix, not - # ScalMat or PDDiagMat (as the `constant` part will have) adaptive = MvNormal(PDMat(Symmetric(dist.Σ))) AdaptiveMvNormal(n, -1, β, σ, dist, adaptive, zeros(n), zeros(n,n)) end +is_symmetric_proposal(::AdaptiveMvNormal) = true + """ adapt!(p::AdaptiveMvNormal, x::AbstractVector) Adaptation for the adaptive multivariate normal mixture proposal as described -in Haario et al. and Roberts & Rosenthal (2009). The code for this routine is -largely taken from `Mamba.jl`. - -# References - -Roberts, Gareth O., and Jeffrey S. Rosenthal. "Examples of adaptive MCMC." -Journal of Computational and Graphical Statistics 18.2 (2009): 349-367. +in Haario et al. (2001) and Roberts & Rosenthal (2009). Will perform an online +estimation of the target covariance matrix and mean. The code for this routine +is largely based on `Mamba.jl`. """ function adapt!(p::AdaptiveMvNormal, x::AbstractVector) p.n += 1 @@ -46,15 +53,13 @@ function adapt!(p::AdaptiveMvNormal, x::AbstractVector) p.EX = f * p.EX + (1.0 - f) * x * x' # compute adapted covariance matrix Σ = (p.σ^2 / (p.d * f)) * (p.EX - p.Ex * p.Ex') - #F = cholesky(Hermitian(Σ), Val{true}(), check=false) F = cholesky(Hermitian(Σ), check=false) if rank(F.L) == p.d - #p.adaptive = MvNormal(PDMat(Symmetric(F.P * F.L))) p.adaptive = MvNormal(PDMat(Σ, F)) end end -function propose(rng::Random.AbstractRNG, p::AdaptiveMvNormal) +function Base.rand(rng::Random.AbstractRNG, p::AdaptiveMvNormal) return if p.n > 2p.d p.β * rand(rng, p.constant) + (1.0 - p.β) * rand(rng, p.adaptive) else @@ -62,8 +67,8 @@ function propose(rng::Random.AbstractRNG, p::AdaptiveMvNormal) end end -function propose(rng::Random.AbstractRNG, p::AdaptiveMvNormal, m::DensityModel) - return propose(rng, p) +function propose(rng::Random.AbstractRNG, proposal::AdaptiveMvNormal, m::DensityModel) + return rand(rng, proposal) end function propose( @@ -73,10 +78,7 @@ function propose( t ) adapt!(proposal, t) - return t + propose(rng, proposal) + return t + rand(rng, proposal) end -# this is always symmetric, so ignore -function q(proposal::AdaptiveMvNormal, t, t_cond) - return 0. -end + diff --git a/test/runtests.jl b/test/runtests.jl index 6b0bd8b..82785b8 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -103,16 +103,16 @@ include("util.jl") end @testset "Adaptive MvNormal mixture ESS" begin - Random.seed!(12) d = 25 - M = randn(d,d) + M = rand(MvNormal(d,1), d) Σ = M*M' m = DensityModel(x -> logpdf(MvNormal(Σ), x)) p = AdaptiveMvNormal(MvNormal(d, 1.)) kwargs = (progress=false, chain_type=Chains) c1 = sample(m, MetropolisHastings(p), 10000; kwargs...) + display(p.adaptive.Σ) c2 = sample(m, RWMH(MvNormal(zeros(d), 1)), 10000; kwargs...) - @test sum(ess(c1).nt.ess .> ess(c2).nt.ess) > 20 + @test sum(ess(c1).nt.ess .> ess(c2).nt.ess) == 25 end @testset "parallel sampling" begin @@ -182,21 +182,22 @@ include("util.jl") @test AdvancedMH.is_symmetric_proposal(p1) # Sample from the posterior with initial parameters. - chain1 = sample(m1, MetropolisHastings(p1), 100000; - chain_type=StructArray, param_names=["x"]) + chain1 = sample(m1, MetropolisHastings(p1), 100000; + progress=false, chain_type=StructArray, param_names=["x"]) @test mean(chain1.x) ≈ mean(d1) atol=0.05 @test std(chain1.x) ≈ std(d1) atol=0.05 end @testset "MALA" begin - + # Set up the sampler. sigma = 1e-1 spl1 = MALA(x -> MvNormal((sigma^2 / 2) .* x, sigma)) # Sample from the posterior with initial parameters. - chain1 = sample(model, spl1, 100000; init_params=ones(2), chain_type=StructArray, param_names=["μ", "σ"]) + chain1 = sample(model, spl1, 100000; progress=false, init_params=ones(2), + chain_type=StructArray, param_names=["μ", "σ"]) @test mean(chain1.μ) ≈ 0.0 atol=0.1 @test mean(chain1.σ) ≈ 1.0 atol=0.1 From 7aa86317b7b160ab1dbbb40d870eaf4184aee1c0 Mon Sep 17 00:00:00 2001 From: Arthur Zwaenepoel Date: Fri, 11 Dec 2020 18:26:09 +0100 Subject: [PATCH 25/31] Update src/adaptivemvnormal.jl Co-authored-by: David Widmann --- src/adaptivemvnormal.jl | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/src/adaptivemvnormal.jl b/src/adaptivemvnormal.jl index bf0372d..19779a2 100644 --- a/src/adaptivemvnormal.jl +++ b/src/adaptivemvnormal.jl @@ -48,9 +48,9 @@ is largely based on `Mamba.jl`. function adapt!(p::AdaptiveMvNormal, x::AbstractVector) p.n += 1 # adapt mean vector and scatter matrix - f = p.n / (p.n + 1.0) - p.Ex = f * p.Ex + (1.0 - f) * x - p.EX = f * p.EX + (1.0 - f) * x * x' + f = p.n / (p.n + 1) + p.Ex = f * p.Ex + (1 - f) * x + p.EX = f * p.EX + (1 - f) * x * x' # compute adapted covariance matrix Σ = (p.σ^2 / (p.d * f)) * (p.EX - p.Ex * p.Ex') F = cholesky(Hermitian(Σ), check=false) @@ -81,4 +81,3 @@ function propose( return t + rand(rng, proposal) end - From 499934999c99307ba27ac36f0c0476b3118d7cdb Mon Sep 17 00:00:00 2001 From: Arthur Zwaenepoel Date: Fri, 11 Dec 2020 18:30:15 +0100 Subject: [PATCH 26/31] Update src/adaptivemvnormal.jl Co-authored-by: David Widmann --- src/adaptivemvnormal.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/adaptivemvnormal.jl b/src/adaptivemvnormal.jl index 19779a2..0508300 100644 --- a/src/adaptivemvnormal.jl +++ b/src/adaptivemvnormal.jl @@ -31,7 +31,7 @@ end function AdaptiveMvNormal(dist::MvNormal; σ=2.38, β=0.05) n = length(dist) - adaptive = MvNormal(PDMat(Symmetric(dist.Σ))) + adaptive = MvNormal(cov(dist)) AdaptiveMvNormal(n, -1, β, σ, dist, adaptive, zeros(n), zeros(n,n)) end @@ -80,4 +80,3 @@ function propose( adapt!(proposal, t) return t + rand(rng, proposal) end - From d4b3f6b705712ad0ec8b5cd9e1b9c9e930416c31 Mon Sep 17 00:00:00 2001 From: Arthur Zwaenepoel Date: Fri, 11 Dec 2020 18:37:15 +0100 Subject: [PATCH 27/31] Update src/adaptivemvnormal.jl Co-authored-by: David Widmann --- src/adaptivemvnormal.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/adaptivemvnormal.jl b/src/adaptivemvnormal.jl index 0508300..aea9f42 100644 --- a/src/adaptivemvnormal.jl +++ b/src/adaptivemvnormal.jl @@ -61,7 +61,7 @@ end function Base.rand(rng::Random.AbstractRNG, p::AdaptiveMvNormal) return if p.n > 2p.d - p.β * rand(rng, p.constant) + (1.0 - p.β) * rand(rng, p.adaptive) + p.β * rand(rng, p.constant) + (1 - p.β) * rand(rng, p.adaptive) else rand(rng, p.constant) end From cb52c7fe331a4f20506062d256293e62f5c54e52 Mon Sep 17 00:00:00 2001 From: Arthur Zwaenepoel Date: Fri, 11 Dec 2020 18:38:59 +0100 Subject: [PATCH 28/31] Update test/Project.toml Co-authored-by: David Widmann --- test/Project.toml | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/test/Project.toml b/test/Project.toml index 5b194a1..8b13789 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -1,4 +1 @@ -[deps] -AdvancedMH = "5b7e9947-ddc0-4b3f-9b55-0d8042f74170" -MCMCChains = "c7f686f2-ff18-58e9-bc7b-31028e88f75d" -StructArrays = "09ab397b-f2b6-538f-b94a-2f83cf4a842a" + From 5a2a175b6db48b15cb4c884e92a5c625f1b11911 Mon Sep 17 00:00:00 2001 From: Arthur Zwaenepoel Date: Fri, 11 Dec 2020 18:43:12 +0100 Subject: [PATCH 29/31] Update src/mh-core.jl Co-authored-by: David Widmann --- src/mh-core.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/mh-core.jl b/src/mh-core.jl index 1f84992..dd32f85 100644 --- a/src/mh-core.jl +++ b/src/mh-core.jl @@ -247,5 +247,4 @@ function AbstractMCMC.step( end end -function accepted!(p::P) where P<:Proposal end - +accepted!(::Proposal) = nothing From b2be967aa1249bddddec5d7c5da7e72dd18628bf Mon Sep 17 00:00:00 2001 From: Arthur Zwaenepoel Date: Fri, 11 Dec 2020 18:43:54 +0100 Subject: [PATCH 30/31] Update src/adaptivemvnormal.jl Co-authored-by: David Widmann --- src/adaptivemvnormal.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/adaptivemvnormal.jl b/src/adaptivemvnormal.jl index aea9f42..de7734d 100644 --- a/src/adaptivemvnormal.jl +++ b/src/adaptivemvnormal.jl @@ -60,7 +60,7 @@ function adapt!(p::AdaptiveMvNormal, x::AbstractVector) end function Base.rand(rng::Random.AbstractRNG, p::AdaptiveMvNormal) - return if p.n > 2p.d + return if p.n > 2 * p.d p.β * rand(rng, p.constant) + (1 - p.β) * rand(rng, p.adaptive) else rand(rng, p.constant) From 518aab10af9bea2377a5b14168d7d5d3a183dbdd Mon Sep 17 00:00:00 2001 From: arzwa Date: Fri, 11 Dec 2020 18:53:53 +0100 Subject: [PATCH 31/31] removed test Project.toml --- test/Project.toml | 1 - 1 file changed, 1 deletion(-) delete mode 100644 test/Project.toml diff --git a/test/Project.toml b/test/Project.toml deleted file mode 100644 index 8b13789..0000000 --- a/test/Project.toml +++ /dev/null @@ -1 +0,0 @@ -