Density Tempered SMC

Given a state space model with a transition vector x, measurement vector y, a vector of parameters for the model θ. The transition density is represented by p(x[t]|x[t-1]) and the measurement density is represented by p(y[t]|x[t]), while the parameters are chosen from a prior density represented by p(θ).

# T = length of vector y
# M = num of θ particles
# N = num of state particles

θ = rand(p(θ_0),M)
θ = Particles(θ)

# initialize X (localized in code to take up less memory)
X = Matrix{Float64}(M,T)

# weight particles by normalizing constant
for m in 1:M
    X[m,:] = bootstrapFilter(N,y,θ)
    θ.logμ[m] = sum(X[m,t].logμ for t in 1:T)
end

ξ = 1.e-12

while ξ ≤ 1.0
    newξ = gridSearch(ξ,θ,B)
    θ = reweight!(θ,(newξ-ξ).*θ.logw)
    
    if θ.ess < B*M
        θ = resample(θ)
        θ = randomwalkMH(θ,y,p(θ))
    end
end

The grid search solves the value of ξ that achieves an effective sample size approximately equal to B*M. This process is sped up by utilizing a bisection method to solve a simple root finding problem. We plug in the function Δ -> ESS(Δ*θ.logw)-B*M and return ξ+Δ. The bisection method code is demonstrated below:

function bisection(f::Function,lower_bound::Float64,upper_bound::Float64,max_iter::Int64=100)
    f_lower_bound = f(lower_bound)
    i = 0

    while upper_bound-lower_bound > 1.e-12
        i += 1
        i != max_iter || error("Max iteration exceeded")

        midpoint = (lower_bound+upper_bound)*0.5
        f_midpoint = f(midpoint)

        if f_midpoint == 0
            break
        elseif f_lower_bound*f_midpoint > 0
            lower_bound = midpoint
            f_lower_bound = f_midpoint
        else
            upper_bound = midpoint
        end
    end

    return upper_bound
end