AdaptiveMCMC

General Usage Performance
,
1

I wrote a Metropolis-Hastings code manually for a problem that samples arrays from a custom distribution that can be found here https://github.com/marouanehanhasse/QComm/blob/main/distributions.jl, before i was using Klara package for older versions of julia, since it was closed i’ve been struggling to find a similar package i tried AdvancedMH but i can’t have outputs since i can’t get the gradient of the target function because it has JuMP, i found AdaptiveMCMC https://juliapackages.com/p/adaptivemcmc and it looks exactly like what i want, but they are giving an example that i didn’t quite understand it for a custom distribution. I want to know how can i use this package if i have a custom distribution from which i can have the proposed array, and a target function.
This is the code where i wrote the manual MCMC:

using Plots
using Random
using SymPy
using CSV
using DataFrames
using Main.QuantumRelay
using Distributions

#calling the github code 
function quantumrelay(alpha,delta,name::String,n)   #n is the number of sources for 1 relay we have 2 sources
   
    chi = fill(sqrt(0.06), n)                  # the parameter chi 
    phi = im*tanh.(chi)
    omega = -log(prod(cosh.(chi)))
    syms, op = qrelay_op(n, phi, alpha, delta)  #operators.jl
    op_a, op_ab, mat, coef = op_mat(op)         #operators.jl

    op_q2 = [syms.apH[1], syms.apV[1], syms.bpH[end], syms.bpV[end]]    # array in operators.jl
    op_q1 = [syms.apH[2:end]..., syms.apV[2:end]..., syms.bpH[1:end-1]..., syms.bpV[1:end-1]...]      #array in operators.jl
    mask_q1 = [op in op_q1 for op in op_a];      #array in operators.jl
    
    mask_q2 = [op in op_q2 for op in op_a];       #array in operators.jl
    qq = [x in syms.apH || x in syms.bpV ? 1 : 0 for x in op_a];      #array in operators.jl
           
    pdet0 = pdet_maker(0.04, 3e-5)      #function utility.jl
    qrs = QRelaySampler(mat, coef, omega, pdet0).         #QuantumRelay.jl  
    targetcache=Dict{Vector{Int}, Float64}()
    target(x::Vector)= log(qrs.prob(qq, x, mask_q1))+ log(qrs.prob(x))      #the log target function of MCMC the function prob is in QuantumRelay.jl 
    Iteration=2^17
    burnin=2^10
    samples=Iteration+burnin
    step=5
    save_iter=burnin:samples
    dist= qrs.psetproposal           #the proposal distribution from the distribution.jl     
    selected=Array[]
    Q = Float64[]
    
    current_x = zeros(8)
    @time for i in 2:samples            #from this line the MCMC algorithm starts
        
        proposed_x= rand(dist(current_x))
        
        prop_proposed= logpdf(dist(current_x), proposed_x)        #the probability from x->y
        
        prop_current= logpdf(dist(proposed_x), current_x)          #the probability from y->x
        
        
        A= min(0,target(proposed_x)+prop_current-target(current_x)-prop_proposed)
        
        
        if log(rand()) < A
            #if i in save_iter
                
                
            push!(selected, proposed_x)
            push!(Q, qrs.prob(qq, proposed_x, mask_q2))
            
            #end
    
            current_x = proposed_x
        else
            current_x= current_x
        end        
    
        
    end
    
    return selected, Q
end

prob= Float64[]
accepted= Array[]
for i = 0:14
    beta = i*pi/14
    name = string(i)
    selected, Q = quantumrelay(pi/4, beta, name,2)
    push!(accepted,length(selected))
    println("beta:", beta)
    push!(prob,Q)
    df=DataFrame(selected=selected)

end

The reason why i want to change it is because using Klara.jl before i was getting better results, and with this one the outcome is not as good as before when i used Klara.jl.

2

Hi,

The main thing you need to do is provide a function of your model’s joint log likelihood. Here is a simple example based on a normal distribution.


using AdaptiveMCMC, Distributions, Random

Random.seed!(845)

function joint_loglike(data, parms)
    μ,σ = parms 
    σ < 0 ? (return -Inf) : nothing
    # log prior of μ
    LL = logpdf(Normal(0, 1), μ)
    # log prior of σ
    LL += logpdf(truncated(Normal(0, 1), 0, Inf), σ)
    # log likelihood of data given μ, σ
    LL += sum(logpdf.(Normal(μ, σ), data))
    return LL 
end

# data for parameter estimation
data = rand(Normal(0, 1), 50)
# number of samples
n = 100_000
L = 2
# initial parameters
init = [0.0,1.0]
output = adaptive_rwm(init, x -> joint_loglike(data, x), div(n,L); L)
samples = output.X'
[mean(samples, dims=1) std(samples, dims=1)]

output

1×4 Matrix{Float64}:
 -0.0154051  0.988926  0.139138  0.101883
2 Likes
3

the distrubution i have i can generate new proposed vector from and i can compute the probability from x->y and from y->x, as for the log target function it’s another function seperated from it.

4

Are you trying to sample from a model with AdaptiveMCMC.jl or are you trying to use a custom sampler with AdaptiveMCMC.jl 's parallel tempering algorithm? If you are trying to do the later, there is an example in a closed issue here.

2 Likes
5

I actually can put anything that can accept my custom proposal distribution that is here https://github.com/marouanehanhasse/QComm/blob/main/distributions.jl, like this proposal function and a log target function which one is more adaptive, because if you take a look at my proposal it returns a type and not a function.

6

in your answer you have your target function join_loglike and your proposed values come from a Normal distribution, what if my proposed values come from a custome proposal distribution such as

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))
7

Can you clarify how you want to use OrthoNNDist? In my example, the proposals are generated by the algorithm in AdaptiveMCMC.jl. The prior distribution of μ is Gaussian, and the distribution of the data is Gaussian.

8

like this distribution generates the proposed values like i said and it can be called for example with this:
x= zeros(8); mat = [0.0 1.0 0.0 0.0 1.0 1.0 0.0 0.0 1.0 0.0 0.0 0.0; 1.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0; 1.0 0.0 1.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 1.0 0.0; 0.0 1.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 1.0 0.0; 0.0 0.0 0.0 1.0 0.0 0.0 0.0 1.0 1.0 0.0 0.0 0.0; 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0; 0.0 0.0 1.0 0.0 0.0 1.0 0.0 0.0 0.0 1.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 1.0 0.0 0.0]
you can call it
OrthoNNDist(x,mat) and you get the following:

OrthoNNDist(
x0: [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
oc: [0.0 1.0 … 0.0 0.0; 1.0 0.0 … 0.0 1.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]
x1s: Any[[0.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0], [1.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0], [0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 0.0], [1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0], [0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0], [0.0, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0, 0.0], [1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0], [0.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0], [0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0]]
prob: 0.08333333333333333
)

if you want to generate an array from x1s which stores all the possible proposal arrays :
rand(OrthoNNDist(x, mat)) and you get one of the arrays stored in x1s:
8-element Vector{Float64}: 1.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0
you can compute pdf and logpdf

9

mainly i want to use this distribution just to generate the proposed arrays and to compute p(x|y) and p(y|x) with y is the proposed value. this distribution used to work fine with Klara.jl

10

I think the best resource for adding a custom sampler to generate proposals can be found here. Aside from that, I am not sure about the details.