-
Notifications
You must be signed in to change notification settings - Fork 0
/
adaptiveMetropolis.jl
82 lines (79 loc) · 2.7 KB
/
adaptiveMetropolis.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
## Makes an adaptive Metropolis kernel proposed by:
## Haario, H., Saksman, E. and Tamminen, J., 2001. An adaptive Metropolis
## algorithm. Bernoulli, 7(2), pp.223-242.
## Essentially provides an adaptive mechanism to make use of the optimal
## scaling results for random walk Metropolis on "Gaussian-like" targets
## described by:
## Roberts, G.O. and Rosenthal, J.S., 2001. Optimal scaling for various
## Metropolis--Hastings algorithms. Statistical Science, 16(4), pp.351-367.
function makeAMKernel(logTargetDensity::F, Σ::SMatrix{d, d, Float64},
updateFrequency::Int64 = 1, ϵ::Float64 = 1.0) where {F<:Function, d}
S::MMatrix{d, d, Float64} = Σ
# A::MMatrix{d, d, Float64} = chol(Symmetric(S))'
A::MMatrix{d, d, Float64} = mychol(S)
scratchv::MVector{d, Float64} = MVector{d, Float64}(undef)
scratchz::MVector{d, Float64} = MVector{d, Float64}(undef)
prevx::MVector{d, Float64} = MVector{d, Float64}(undef)
ldprevx = Ref(-Inf)
accepts = Ref(0)
calls = Ref(0)
covEstimate::MMatrix{d, d, Float64} = Σ
meanEstimate::MVector{d, Float64} = zeros(MVector{d, Float64})
@inline function retuneSigma()
if covEstimate[1,1] == 0.0
S .= Σ * ϵ / calls.x
else
S .= 5.6644 / d * covEstimate
end
try
# A .= chol(Symmetric(S))'
A .= mychol(S)
## quick fix as chol on SMatrix doesn't throw not pos def exceptions
any(isnan, A) && throw(DomainError())
catch e
S .= Σ * ϵ / calls.x
# A .= chol(Symmetric(S))'
A .= mychol(S)
end
end
@inline function P(x::SVector{d, Float64})
calls.x += 1
randn!(scratchv)
mul!(scratchz, A, scratchv)
z::SVector{d, Float64} = scratchz + x
# scratchv .= A * scratchv
# z::SVector{d, Float64} = x + scratchv
if x == prevx
lpi_x = ldprevx.x
else
lpi_x = logTargetDensity(x)
prevx .= x
end
lpi_z = logTargetDensity(z)
if -randexp() < lpi_z - lpi_x
prevx .= z
ldprevx.x = lpi_z
accepts.x += 1
rval = z
else
rval = x
end
t::Int64 = calls.x
covEstimate .= (t-1)/t * (covEstimate +
(rval - meanEstimate) * (rval - meanEstimate)' / t)
meanEstimate .= (t-1)/t.*meanEstimate + rval/t
mod(t, updateFrequency) == 0 && retuneSigma()
return rval
end
@inline function P(s::Symbol)
s == :acceptanceRate && return accepts.x / calls.x
s == :meanEstimate && return meanEstimate
s == :covEstimate && return covEstimate * calls.x /(calls.x-1)
end
return P
end
function makeAMKernel(logTargetDensity::F, d::Int64, updateFrequency::Int64 = 1,
ϵ::Float64 = 1.0) where F<:Function
Id = SMatrix{d, d, Float64}(Matrix(1.0I, d, d))
return makeAMKernel(logTargetDensity::F, Id, updateFrequency, ϵ)
end