How to call a customized proposal function in AdaptiveMCMC.jl

Viewed 32

I've been trying to use a customized proposal distribution that generates the proposed arrays so we can use them and test or sample them in a Metropolis-Hastings algorithm, with a log_target function, i wrote the metropolis-hastings code manually and it works fine, tho it doesn't give the same satisfactory results as using Klara.jl which is a closed package now. The proposal distribution is written like this

using Distributions
nonneg(v) = all(v.>=0) ? true : false

struct OrthoNNDist <: DiscreteMultivariateDistribution
    x0::Vector{Float64}
    oc::Array{Float64,2}
    x1s::Array
    prob::Float64
    #return a new uniform distribution with all vectors in x1s orthogonal to oc
    function OrthoNNDist(x0::Vector{Float64}, oc::Array{Float64,2})
        x1s = []
        for i = 1:size(oc)[2]
            x1 = x0 + oc[:, i]
            if nonneg(x1)
                push!(x1s, x1)
            end
            x1 = x0 - oc[:, i]
            if nonneg(x1)
                push!(x1s, x1)
            end
        end
        new(x0, oc, x1s, 1.0/length(x1s))
    end
end

Base.length(d::OrthoNNDist) = length(d.x0)

Distributions.rand(d::OrthoNNDist) = rand(d.x1s)

Distributions.pdf(d::OrthoNNDist, x::Vector) = x in d.x1s ? d.prob : 0.0
Distributions.pdf(d::OrthoNNDist) = fill(d.prob, size(d.x1s))
Distributions.logpdf(d::OrthoNNDist, x::Vector) = log(pdf(d, x))

you can test it using for example x0=[1.0, 1.0, 0.0,0.0, 0.0, 0.0, 1.0, 1.0] and mat = [0 1 0 0 1 1 0 0 1 0 0 0; 1 0 0 0 1 0 0 0 0 0 0 1; 1 0 1 0 0 0 0 1 0 0 1 0; 0 1 0 0 0 0 1 0 0 0 1 0; 0 0 0 1 0 0 0 1 1 0 0 0; 0 0 0 1 0 0 0 0 0 0 0 1; 0 0 1 0 0 1 0 0 0 1 0 0; 0 0 0 0 0 0 1 0 0 1 0 0]. we can fix the mat and let x0 to be the variable by writing for example proposal(x::Vector)= OrthoNNDist(x,mat). If we have a log_target function in general called logtarget(x::Vector) and using this proposal distribution above and i want to use it in AdaptiveMCMC package or any other package that can be used in this case with a minimum example, i tried AdvancedMH but a part of my target function uses JuMP hence i can't get the gradient of the target, i've read the Mamba documentation but i couldn't understand it exactly i would with a minimum example in this case of a customized proposal function and a target function, AdaptiveMCMC looks more simple but i keep getting some MethodError regarding the distribution it's supposed to be a function but mine it's not as you can see in the code above. i can provide here an example of the code using Klara with comments.

function qrelay(alpha, delta, name)
    n = 2
    chi = fill(sqrt(0.06), n)
    phi = im * tanh(chi)
    omega = 1.0 / prod(cosh(chi))^2
    syms, op = qrelay_op(n, phi, alpha, delta) #it gives an array
    op_a, op_ab, mat, coef = op_mat(op)  #array, array, matrice, array of coefficients

    op_q2 = [syms.apH[1], syms.apV[1], syms.bpH[end], syms.bpV[end]] #array
    op_q1 = [syms.apH[2:end]..., syms.apV[2:end]..., syms.bpH[1:end-1]...,  syms.bpV[1:end-1]...] #array
    mask_q1 = [op in op_q1 for op in op_a]; #array
    mask_q2 = [op in op_q2 for op in op_a];  #array
    qq = [x in syms.apH || x in syms.bpV ? 1 : 0 for x in op_a] #array
    
    pdet0 = pdet_maker(0.04, 1e-5) #it gives a probability
    qrs = QSampler(mat, coef, omega, pdet0) #calling a module QSampler

    targetcache = Dict{Vector{Int}, Float64}()
    function plogtarget(na::Vector{Int})
        get!(targetcache, na) do
            log(qrs.prob(qq, na, mask_q1) * qrs.prob(na))
        end
    end
#     plogtarget(na::Vector{Int}) = log(qrs.prob(qq, na, mask_q1) * qrs.prob(na))
    p = BasicDiscMuvParameter(:p, logtarget=plogtarget)
    model = likelihood_model([p], isindexed=false)
    sampler = MH(qrs.psetproposal, symmetric=false)    #this is where the proposal function is called qrs.psetproposal(x::Vector)= qrs.OthoNNDist(x, mat)
    mcrange = BasicMCRange(nsteps=2^20 + 2^10, burnin=2^10, thinning=2^5)
    v0 = Dict(:p=>zeros(qq))
    outopts = Dict{Symbol, Any}(
        :monitor=>[:value, :logtarget],
        :diagnostics=>[:accept],
        :destination=>:iostream,
        :filepath=>"$dataname/"*name
    )
    job = BasicMCJob(model, sampler, mcrange, v0, outopts=outopts)

    funcQ(v) = qrs.prob(qq, v, mask_q2)  #it returns a probability
    
    return qrs, job, funcQ
end

This piece of code won't work of course because Klara is closed but i just put it here to give you a clearer view of how the code was working before.

0 Answers
Related