This note is mainly based on Chapter 12 of Liu (2008).

where MTM refers to

and its simplified version **orientational bias Monte Carlo (OBMC)** when $w(\mathbf x, \mathbf y)=\pi(\mathbf x)$:

whose contour plot is

For the general MH procedure:

function propose(x; a = 4)θ = rand() * 2πr = rand() * 4return [x[1] + r*cos(θ), x[2] + r*sin(θ)]endfunction iMH(M = 200000)# start pointsx = [rand() - 0.5, rand() - 0.5]X = ones(M, 2)num = 0for m = 1:Mxs = propose(x)ρ = targetpdf(xs) / targetpdf(x)if rand() < ρx = xsnum += 1endX[m,:] = xendreturn X, num/Mend

For CGMC, there is a long way to go. Firstly, find out the derivative of target pdf and find the anchor point:

# derivative of targetpdffunction dm(x; step = 0.05, nmaxstep = 100)term1 = w[1] / 2π * exp(-0.5*sum(x.^2))term2 = w[2] / (2π * 0.19) * exp( -0.5 * transpose(x.+6) * invΣ2 * (x.+6) )term3 = w[3] / (2π * 0.19) * exp( -0.5 * transpose(x.-4) * invΣ3 * (x.-4) )term = term1 .+ term2 * invΣ2 + term3 * invΣ3# negative of gradient direction (remove the minus symbol already)dx1 = term[1,1]*x[1] + term[1,2]*x[2]dx2 = term[2,1]*x[1] + term[2,2]*x[2]dx = [dx1, dx2]# find local maximumbst = targetpdf(x)i = 0while i < nmaxstepi += 1cur = targetpdf(x .- step * dx)if cur > bstbst = curx = x .- step * dxelsebreakendendreturn xend

Then implement MTM algorithm for sampling radius $r$:

function MTM(anchor, direction;M = 100, k = 5, σ = 10)# use the simplified version: OBMC algorithm# initial xx = 2*rand()-1num = 0for m = 1:My = ones(k)wy = ones(k)for i = 1:ky[i] = rand(Normal(x, σ))wy[i] = f(y[i], anchor, direction)end#wy = wy ./ sum(wy)yl = sample(y, pweights(wy))# draw reference pointsxprime = ones(k)xprime[k] = xwxprime = ones(k)wxprime[k] = f(x, anchor, direction)for i = 1:k-1xprime[i] = rand(Normal(yl, σ))wxprime[i] = f(xprime[i], anchor, direction)end# accept or notρ = sum(wy) / sum(wxprime)if rand() < ρx = ylnum += 1endendreturn x, num/Mend

Finally, we can combine these sub-procedures:

function cgmc(M = 200000)# start points (each col is a stream)x = rand(2, 2) .- 0.5X = ones(M, 4)for m = 1:M# randomly choose xaxaIdx = sample([1,2])xa = x[:,xaIdx]# find gradient of pdf at xa (or conjugate gradient of log pdf at xa)y = dm(xa)xcIdx = 3 - xaIdxxc = x[:,xcIdx]e = y - xce = e / norm(e)# sample rr, = MTM(y, e)# update xcx[:,xcIdx] = y .+ r*eX[m, :] .= vec(x)endreturn Xend

Refer to bimodal-example.jl for the complete source code. Run the program, I can reproduce Fig 11.1 successfully,