1+ struct Ensemble{D} <: MHSampler
2+ n_walkers:: Int
3+ proposal:: D
4+ end
5+
6+ # Define the first step! function, which is called at the
7+ # beginning of sampling. Return the initial parameter used
8+ # to define the sampler.
9+ function AbstractMCMC. step! (
10+ rng:: Random.AbstractRNG ,
11+ model:: DensityModel ,
12+ spl:: Ensemble ,
13+ N:: Integer ,
14+ :: Nothing ;
15+ init_params = nothing ,
16+ kwargs... ,
17+ )
18+ if init_params === nothing
19+ return propose (rng, spl, model)
20+ else
21+ return Transition (model, init_params)
22+ end
23+ end
24+
25+ # Define the other step functions. Returns a Transition containing
26+ # either a new proposal (if accepted) or the previous proposal
27+ # (if not accepted).
28+ function AbstractMCMC. step! (
29+ rng:: Random.AbstractRNG ,
30+ model:: DensityModel ,
31+ spl:: Ensemble ,
32+ :: Integer ,
33+ params_prev;
34+ kwargs... ,
35+ )
36+ # Generate a new proposal. Accept/reject happens at proposal level.
37+ return propose (rng, spl, model, params_prev)
38+ end
39+
40+ #
41+ # Initial proposal
42+ #
43+ function propose (rng:: Random.AbstractRNG , spl:: Ensemble , model:: DensityModel )
44+ # Make the first proposal with a static draw from the prior.
45+ static_prop = StaticProposal (spl. proposal. proposal)
46+ mh_spl = MetropolisHastings (static_prop)
47+ return [propose (rng, mh_spl, model) for _ in 1 : spl. n_walkers]
48+ end
49+
50+ #
51+ # Every other proposal
52+ #
53+ function propose (rng:: Random.AbstractRNG , spl:: Ensemble , model:: DensityModel , walkers:: Vector{W} ) where {W<: Transition }
54+ new_walkers = similar (walkers)
55+
56+ others = 1 : (spl. n_walkers - 1 )
57+ for i in 1 : spl. n_walkers
58+ walker = walkers[i]
59+ idx = mod1 (i + rand (rng, others), spl. n_walkers)
60+ other_walker = walkers[idx]
61+ new_walkers[i] = move (rng, spl, model, walker, other_walker)
62+ end
63+
64+ return new_walkers
65+ end
66+
67+
68+ # ####################################
69+ # Basic stretch move implementation #
70+ # ####################################
71+ struct StretchProposal{P, F<: AbstractFloat } <: Proposal{P}
72+ proposal :: P
73+ stretch_length:: F
74+ end
75+
76+ StretchProposal (p) = StretchProposal (p, 2.0 )
77+
78+ function move (
79+ rng:: Random.AbstractRNG ,
80+ spl:: Ensemble{<:StretchProposal} ,
81+ model:: DensityModel ,
82+ walker:: Transition ,
83+ other_walker:: Transition ,
84+ )
85+ # Calculate intermediate values
86+ proposal = spl. proposal
87+ n = length (walker. params)
88+ a = proposal. stretch_length
89+ z = ((a - 1 ) * rand (rng) + 1 )^ 2 / a
90+ alphamult = (n - 1 ) * log (z)
91+
92+ # Make new parameters
93+ y = @. walker. params + z * (other_walker. params - walker. params)
94+
95+ # Construct a new walker
96+ new_walker = Transition (model, y)
97+
98+ # Calculate accept/reject value.
99+ alpha = alphamult + new_walker. lp - walker. lp
100+
101+ if - Random. randexp (rng) <= alpha
102+ return new_walker
103+ else
104+ return walker
105+ end
106+ end
107+
108+ # ########################
109+ # Elliptical slice step #
110+ # #########################
111+
112+ # struct EllipticalSlice{E} <: ProposalStyle
113+ # ellipse::E
114+ # end
115+
116+ # function move(
117+ # # spl::Ensemble,
118+ # spl::Ensemble{Proposal{T,P}},
119+ # model::DensityModel,
120+ # walker::Transition,
121+ # other_walker::Transition,
122+ # ) where {T<:EllipticalSlice,P}
123+ # # Calculate intermediate values
124+ # proposal = spl.proposal
125+ # n = length(walker.params)
126+ # nu = rand(proposal.type.ellipse)
127+
128+ # u = rand()
129+ # y = walker.lp - Random.randexp()
130+
131+ # theta = 2 * π * rand()
132+
133+ # theta_min = theta - 2.0*π
134+ # theta_max = theta
135+
136+ # f = walker.params
137+ # while true
138+ # stheta, ctheta = sincos(theta)
139+
140+ # f_prime = f .* ctheta + nu .* stheta
141+
142+ # new_walker = Transition(model, f_prime)
143+
144+ # if new_walker.lp > y
145+ # return new_walker
146+ # else
147+ # if theta < 0
148+ # theta_min = theta
149+ # else
150+ # theta_max = theta
151+ # end
152+
153+ # theta = theta_min + (theta_max - theta_min) * rand()
154+ # end
155+ # end
156+ # end
157+
158+ # ####################
159+ # Slice and stretch #
160+ # ####################
161+ # struct EllipticalSliceStretch{E, S<:Stretch} <: ProposalStyle
162+ # ellipse::E
163+ # stretch::S
164+ # end
165+
166+ # EllipticalSliceStretch(e) = EllipticalSliceStretch(e, Stretch(2.0))
167+
168+ # function move(
169+ # # spl::Ensemble,
170+ # spl::Ensemble{Proposal{T,P}},
171+ # model::DensityModel,
172+ # walker::Transition,
173+ # other_walker::Transition,
174+ # ) where {T<:EllipticalSliceStretch,P}
175+ # # Calculate intermediate values
176+ # proposal = spl.proposal
177+ # n = length(walker.params)
178+ # nu = rand(proposal.type.ellipse)
179+
180+ # # Calculate stretch step first
181+ # subspl = Ensemble(spl.n_walkers, Proposal(proposal.type.stretch, proposal.proposal))
182+ # walker = move(subspl, model, walker, other_walker)
183+
184+ # u = rand()
185+ # y = walker.lp - Random.randexp()
186+
187+ # theta = 2 * π * rand()
188+
189+ # theta_min = theta - 2.0*π
190+ # theta_max = theta
191+
192+ # f = walker.params
193+
194+ # i = 0
195+ # while true
196+ # i += 1
197+
198+ # stheta, ctheta = sincos(theta)
199+
200+ # f_prime = f .* ctheta + nu .* stheta
201+
202+ # new_walker = Transition(model, f_prime)
203+
204+ # # @info "Slice step" i f f_prime y new_walker.lp theta theta_max theta_min
205+
206+ # if new_walker.lp > y
207+ # return new_walker
208+ # else
209+ # if theta < 0
210+ # theta_min = theta
211+ # else
212+ # theta_max = theta
213+ # end
214+
215+ # theta = theta_min + (theta_max - theta_min) * rand()
216+ # end
217+ # end
218+ # end
0 commit comments