blaudichte <- function(x) {
.3333*dnorm(x[1],2,.5)*dnorm(x[2],3)+
.3333*dnorm(x[1],5,1)*dnorm(x[2],1,1/3)+
.3333*dnorm(x[1],3,0.5)*dnorm(x[2],1.5,.5)
}

a <- rep(1:80/10,80)
b <- rep(1:80/10,rep(80,80))
contour(1:80/10,1:80/10,matrix(apply(cbind(a,b),1,blaudichte),80),col="blue",nlevels=20)

mcmcsample<-function(n,b=0,s=1) {
x <- c(7,7)
a<-c()
b<-c()
j<-0
for (i in 1:n) {
y <- x+rnorm(2)/20
if (i%%s==0) {
j<-j+1
a[j]<-x[1]
b[j]<-x[2]
}
if (blaudichte(y)/blaudichte(x)>runif(1)) {
lines(c(x[1],y[1]),c(x[2],y[2]),col='red') 
x<-y
}
}
return(list(a,b))
}

a <- rep(1:80/10,80)
b <- rep(1:80/10,rep(80,80))
contour(1:80/10,1:80/10,matrix(apply(cbind(a,b),1,blaudichte),80),col="blue",nlevels=20)

mcmcsample(30000) -> erg 


contour(1:80/10,1:80/10,matrix(apply(cbind(a,b),1,blaudichte),80),col="blue",nlevels=20)
points(erg[[1]],erg[[2]],cex=0.2)


contour(1:80/10,1:80/10,matrix(apply(cbind(a,b),1,blaudichte),80),col="blue",nlevels=20)
points(erg[[1]][2000:30000],erg[[2]][2000:30000],cex=0.2)


contour(1:80/10,1:80/10,matrix(apply(cbind(a,b),1,blaudichte),80),col="blue",nlevels=20)
points(erg[[1]][10:150*200],erg[[2]][10:150*200],cex=0.2)