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.