-
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 Metropolis algorithm #57
base: master
Are you sure you want to change the base?
Changes from 5 commits
71f1072
3a23129
c40d4db
5afb87e
8cb275d
46f77ca
970b0f5
c800aaf
eb91b36
77bd3e5
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,56 @@ | ||
# Simple Adaptive Metropolis Proposal | ||
# The proposal at each step is equal to a scalar multiple | ||
# of the empirical posterior covariance plus a fixed, small covariance | ||
# matrix epsilon which is also used for initial exploration. | ||
# | ||
# Reference: | ||
# H. Haario, E. Saksman, and J. Tamminen, "An adaptive Metropolis algorithm", | ||
# Bernoulli 7(2): 223-242 (2001) | ||
mutable struct AMProposal{FT <: Real, CT <: AbstractMatrix{FT}, | ||
MNT <: AbstractMvNormal} <: Proposal{MNT} | ||
epsilon::CT | ||
scalefactor::FT | ||
proposal::MNT | ||
samplemean::Vector{FT} | ||
samplesqmean::Matrix{FT} | ||
N::Int | ||
end | ||
|
||
function AMProposal(epsilon::AbstractMatrix{FT}, | ||
scalefactor=FT(2.38^2 / size(epsilon, 1))) where {FT <: Real} | ||
kaandocal marked this conversation as resolved.
Show resolved
Hide resolved
|
||
sym = PDMat(Symmetric(epsilon)) | ||
proposal = MvNormal(zeros(size(sym, 1)), sym) | ||
AMProposal(sym, scalefactor, proposal, mean(proposal), zeros(size(sym)...), 0) | ||
end | ||
|
||
logratio_proposal_density(::AMProposal, params_prev, params) = 0 | ||
|
||
# When the proposal is initialised the empirical posterior covariance is zero | ||
function trackstep(proposal::AMProposal, params) | ||
proposal.samplemean .= params | ||
proposal.samplesqmean .= params * params' | ||
proposal.proposal = MvNormal(zeros(size(proposal.samplemean)), proposal.epsilon) | ||
proposal.N = 1 | ||
end | ||
|
||
# Recompute the empirical posterior covariance matrix | ||
function trackstep(proposal::AMProposal, params, ::Union{Val{true}, Val{false}}) | ||
proposal.samplemean .= (proposal.samplemean .* proposal.N .+ params) ./ (proposal.N + 1) | ||
|
||
proposal.samplesqmean .= (proposal.samplesqmean .* proposal.N + params * params') ./ | ||
(proposal.N + 1) | ||
|
||
samplecov = proposal.samplesqmean .- proposal.samplemean * proposal.samplemean' | ||
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. One should be aware that this is the "naive" algorithm for computing the variance which can lead to catastrophic cancellation (see e.g. https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Na%C3%AFve_algorithm). Maybe it would be better to use Welford's algorithm (https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Welford's_online_algorithm)? 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. I will have a look at this later. In combination with creating a new 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. I think the Welford implementation in AdvancedHMC is pretty good. I think it's an adaptation of the Stan version. 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. Added! |
||
prop_cov = proposal.scalefactor .* samplecov | ||
proposal.proposal = MvNormal(mean(proposal.proposal), prop_cov .+ proposal.epsilon) | ||
proposal.N += 1 | ||
end | ||
|
||
function propose(rng::Random.AbstractRNG, p::AMProposal, ::DensityModel) | ||
return rand(rng, p.proposal) | ||
end | ||
|
||
function propose(rng::Random.AbstractRNG, p::AMProposal, ::DensityModel, t) | ||
return t + rand(rng, p.proposal) | ||
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.
Can you add a
[compat]
entry for PDMats?