## TARGET N(0,1)
## Gaussian kernel
Niter <- 10000
x <- rep(0,Niter)
x[1] <- rnorm(1)
sigma2 <- 1
for (i in 2:Niter)
{
xtilde <- rnorm(1,x[i-1],sqrt(sigma2))
rho <- dnorm(xtilde)/dnorm(x[i-1])
if (runif(1)<=rho) x[i] <- xtilde else x[i] <- x[i-1]
}
plot(x,type="l",lwd=2,col="red")
plot(acf(x,lag=200))
hist(x,prob=TRUE,col="tomato")
curve(dnorm,add=TRUE,col="blue",lwd=2)
indi <- seq(1,Niter,by=20)
hist(x[indi],prob=TRUE,nclass=30,col="tomato")
curve(dnorm,add=TRUE,col="blue",lwd=2)
## TARGET gamma(2,1)
## Gaussian kernel
Niter <- 10000
x <- rep(0,Niter)
x[1] <- 20
sigma2 <- 8
lulu <- 0
for (i in 2:Niter)
{
xtilde <- rnorm(1,x[i-1],sqrt(sigma2))
if (xtilde>0) rho <- dgamma(xtilde,2,1)/dgamma(x[i-1],2,1)
else rho <- 0
if (runif(1)<=rho) { x[i] <- xtilde ; lulu <- lulu+1 } else x[i] <- x[i-1]
}
plot(x[-(1:500)],type="l",lwd=2,col="red")
plot(acf(x[-(1:500)],lag=200))
hist(x[-(1:500)],prob=TRUE,col="tomato")
curve(dgamma(x,2,1),add=TRUE,col="blue",lwd=2)
## TARGET Gaussian mixture
## Gaussian kernel
target <- function(x) {0.3*dnorm(x,1,1)+
0.7*dnorm(x,12,1)}
curve(target,-5,15,col="blue",lwd=2)
Niter <- 1000000
x <- rep(0,Niter)
x[1] <- 10
sigma2 <- 100
lulu <- 0
for (i in 2:Niter)
{
xtilde <- rnorm(1,x[i-1],sqrt(sigma2))
rho <- target(xtilde)/target(x[i-1])
if (runif(1)<=rho) { x[i] <- xtilde ; lulu <- lulu+1 } else x[i] <- x[i-1]
}
plot(x[10^4:(2*10^4)],type="l",lwd=2,col="red")
plot(acf(x[10^4:(2*10^4)],lag=200))
hist(x[-(1:500)],prob=TRUE,col="tomato")
curve(target,add=TRUE,col="blue",lwd=2)
## TARGET Gaussian mixture
## Independent Gaussian kernel
target <- function(x) {0.3*dnorm(x,1,1)+
0.7*dnorm(x,12,1)}
curve(target,-5,15,col="blue",lwd=2)
Niter <- 1000000
x <- rep(0,Niter)
x[1] <- 10
sigma2 <- 100
lulu <- 0
for (i in 2:Niter)
{
xtilde <- rnorm(1,6.5,sqrt(sigma2))
rho <- (target(xtilde)/dnorm(xtilde,6.5,sqrt(sigma2)))/
(target(x[i-1])/dnorm(x[i-1],6.5,sqrt(sigma2)))
if (runif(1)<=rho) { x[i] <- xtilde ; lulu <- lulu+1 } else x[i] <- x[i-1]
}
plot(x[10^4:(2*10^4)],type="l",lwd=2,col="red")
plot(acf(x[10^4:(2*10^4)],lag=200))
hist(x[10^4:(2*10^4)],prob=TRUE,col="tomato")
curve(target,add=TRUE,col="blue",lwd=2)