We can write the following Julia code which use uniform distribution U[−δ,δ] as g.
# random walk Metropolis-Hastingsfunctionrmh(T, delta, f::Function, initval =5) x =ones(T+1) x[1] = initvalfor t =1:T# generate Yt epsi =rand() *2* delta - delta Yt = epsi + x[t]# accept or not u =rand() r =f(Yt)/f(x[t])if r >=1 x[t+1] = Ytelseif u <= r x[t+1] = Ytelse x[t+1] = x[t]endendendreturn(x)end
Then apply this algorithm to the normal distribution:
# density function of N(0, 1) without normalizationfunctiondnorm(x)return(exp(-0.5* x^2))end# exampledeltalist = [0.1,0.5,1]N =15000# results of mean and variancemuvar =zeros(3,2)using Statisticsps = []for i =1:3 x =rmh(N, deltalist[i], dnorm)push!(ps,plot(x, legend=:none)) muvar[i,1] =mean(x) muvar[i,2] =var(x)end# plotplot(ps[1], ps[2], ps[3], layout=(3,1))savefig("res_rmh.png")