Burn-in in Metropolis-Hastings

Specific Domains Probabilistic programming
1

i have written a metropolis-hastings code for myself but it doesn’t quite give me the smoothest result that i want so i want to put the burn-in in it, maybe finding a good starting point will increase the efficiency, the code sample some vectors that and then i compute a probability of having a specific vector from the results of the code then compute the mean of these probability then plot it, so here’s my code and i’ll appreciate any suggestion to increase the efficiency of the code

using Distributions
using LinearAlgebra
using Plots
using Random
using QuantumRelay      # A customized package
using SymPy
using CSV
using DataFrames

function relay_basis(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 = 1.0/prod(cosh.(chi))^2
    syms, op = qrelay_op(n, phi, alpha, delta)
    op_a, op_ab, mat, coef = op_mat(op)

    op_q2 = [syms.apH[1], syms.apV[1], syms.bpH[end], syms.bpV[end]]
    op_q1 = [syms.apH[2:end]..., syms.apV[2:end]..., syms.bpH[1:end-1]..., syms.bpV[1:end-1]...]
    mask_q1 = [op in op_q1 for op in op_a];
    
     mask_q2 = [op in op_q2 for op in op_a];
     qq = [x in syms.apH || x in syms.bpV ? 1 : 0 for x in op_a]
           
    pdet0 = pdet_maker(0.04, 1e-5)
    qrs = QRelaySampler(mat, coef, omega, pdet0)
    targetcache=Dict{Vector{Int}, Float64}()
    logtarget(x::Vector)= log(qrs.prob(qq, x, mask_q1) * qrs.prob(x))      #the log target function of MCMC
    Iteration=10000
    burnin=500
    samples=Iteration+burnin
    dist= qrs.psetproposal           #the proposal distribution, it takes a vector as an argument
    
    selected=[]
    Q = []
    
    current_x = [0, 0, 0, 0, 0, 0, 0, 0]      #the initial point
    @time for i in 2:Iteration            #from this line the MCMC algorithm starts
        
        proposed_x= rand(dist(current_x))
        a=dist(current_x).prob/dist(proposed_x).prob
        C= min(1,a*exp(target(proposed_x)-target(current_x)))
        
        
        if rand() < C
            push!(selected, proposed_x)
            push!(Q, qrs.prob(qq, proposed_x, mask_q2)) #computing the probability of getting qq knowing proposed_x
            current_x = proposed_x
        end
        
    end
    
    return selected, Q
end

i want to know how can i perform the burn-in inside this function?

2

I like to define a range of samples to actually collect:

save_iter = burnin:step:Iterations

and then just check

   if rand() < C
       if i in save_iter
            push!(selected, proposed_x)
            push!(Q, qrs.prob(qq, proposed_x, mask_q2)) #computing the probability of getting qq knowing proposed_x
       end
       current_x = proposed_x
   end

That way you get burnin and subsampling in one interface

3 Likes
3

I like to plot the trace before doing the burning, because if your starting point is far from the main mode of the distribution it can take a while to converge, so I’d rather do it outside of M-H.

1 Like
4

i didn’t quite understand what you said, how can i plot the trace and what do you mean by do it outside MH

5

Just plot the variable you are sampling (selected in your code) in function of the iterations. Here for example I sample a Normal distribution centered in (0,0) and an initial point at (10,10), as you can see it takes about 200 iterations for the sampler to converge towards the mode of the distribution, so I would burn about 200 iterations (e.g. burn_selected = selected[200:step:end]). Without plotting it’s hard to know if you are burning too little or too much.

tmp
tmp945×787 36.9 KB

(I’m just plotting the first dimension)

3 Likes
6

it is a nice idea except in my code there is a customized distribution that i’m not sure where does it converge

7

thanks for the idea, perhaps i want to know also is there any way that i can perform thinning in the code?

8

What is thinning?

1 Like
9

Thinning = subsampling here. @marouane1994 Just set step in my code above to something larger than 1 and you start to skip samples in push!.

2 Likes
10

thinning means for example that every third accepted value will be removed if thinning=3

1 Like