# Metropolis Adjusted Langevin Algorithm (MALA)

Haven’t dumped much code here in a while. Here’s a Julia implementation of MALA with an arbitrary preconditioning matrix M. Potentially I might use this in the future.

## Generic Julia Implementation

Arguments are a function to evaluate the logdensity, function to evaluate the gradient, a step size h, a preconditioning matrix M, number of iterations and an initial parameter value.

function mala(logdensity,gradient,h,M,niter,θinit)
function gradientStep(θ,t)
θ-t*M*gradient(θ)
end
θtrace=Array{Float64}(length(θinit),niter)
θ=θinit
θtrace[:,1]=θinit
for i=2:niter
θold=θ
θ=rand(MvNormal(gradientStep(θ,0.5*h),h*M))
d=logdensity(θ) - logdensity(θold) + logpdf(MvNormal(gradientStep(θ,0.5*h),h*M),θold) - logpdf(MvNormal(gradientStep(θold,0.5*h),h*M),θ)
if(!(log(rand(Uniform(0,1)))<d))
θ=θold
end
θtrace[:,i]=θ
end
θtrace
end


Don’t forget to use the Distributions package of Julia.

## Bivariate Normal Example

Consider a Bivariate Normal as an example

ρ²=0.8
Σ=[1 ρ²;ρ² 1]

function logdensity(θ)
logpdf(MvNormal(Σ),θ)
end

function gradient(θ)
Σ\θ
end


and now the code to generate the plots and results

niter=1000
h=1/eigs(inv(Σ),nev=1)[1][1]
@time draws=mala(logdensity,gradient,h,eye(2),niter,[5,50]);
sum(map(t -> draws[:,t]!=draws[:,t-1],2:niter))/(niter-1)
@time pdraws=mala(logdensity,gradient,h,Σ,niter,[5,50]);
sum(map(t -> pdraws[:,t]!=pdraws[:,t-1],2:niter))/(niter-1)

function logdensity2d(x,y)
logdensity([x,y])
end
x = -30:0.1:30
y = -30:0.1:50
X = repmat(x',length(y),1)
Y = repmat(y,1,length(x))
Z = map(logdensity2d,Y,X)
p1 = contour(x,y,Z,200)
plot(vec(draws[1,:]),vec(draws[2,:]))
plot(vec(pdraws[1,:]),vec(pdraws[2,:]))


pdraws uses the covariance matrix Σ as the preconditioning matrix, whereas the first uses an identity matrix, resulting in the original MALA algorithm. The traceplot of draws from MALA and preconditioned MALA are shown in blue and green respectively…

Traceplots for MALA and preconditioned MALA

## Effective Sample Sizes in R

We can use the julia “RCall” package to switch over to R and use the coda library to evaluate the minimum effective sample size for both of these MCMC algorithms.

julia> library(coda)

R> library(coda)

R> min(effectiveSize($(draws'))) [1] 22.02418 R> min(effectiveSize($(pdraws')))
[1] 50.85163


I didn’t tune the step size h in this example at all (you should).