# Burn-in in Metropolis-Hastings

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, syms.apV, 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?

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

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

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

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.

(I’m just plotting the first dimension)

3 Likes

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

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

What is thinning?

1 Like

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

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

1 Like