options(repr.plot.width=5, repr.plot.height=4)
library(ggplot2)
library(bayesplot)

targetD <- function(x) {
    return(0.7 * dnorm(x, 0, 1) + 0.3 * dnorm(x, 2, .5))
    }

x <- seq(-5, 5, .01)
plot(x, targetD(x), type='l')

Metropolis-Hastings

stepsize <- 1
proposal <- function(x) {
    return(rnorm(1, x, stepsize))
}

initial <- 0
NSteps <- 10000

mcmc <- function(targetD, initial=0, proposalD=proposal, N=NSteps) {
    chain <- initial
    current <- initial
    for (i in 1:N) {
        new <- proposalD(current)
        if (runif(1) <= min(1, targetD(new)/targetD(current))) {
            current <- new
        }
        chain <- c(chain, current)
    }
    return(chain)
}

chain1 <- mcmc(targetD, initial=-10)
chain2 <- mcmc(targetD, initial=10)

chains <- array(dim = c(length(chain1), 2, 1), 
                dimnames=list(iterations=NULL, chains=c('chain1', 'chain2'), parameters='x'))

chains[,1,1] <- chain1
chains[,2,1] <- chain2

color_scheme_set("red")

mcmc_trace(chains, n_warmup = 1000)

mcmc_areas(chains[-(1:100),,], prob=0.95)

mcmc_acf(chains[-(1:100), , ])

library(coda)

Attaching package: ‘coda’

The following object is masked _by_ ‘.GlobalEnv’:

    mcmc

s1 <- coda::mcmc(data=chain1[1:5000])
s2 <- coda::mcmc(data=chain2[1:5000])

gelman.diag(mcmc.list(list(s1, s2)), autoburnin = T)

Potential scale reduction factors:

     Point est. Upper C.I.
[1,]          1          1