prior_mean <- 0 prior_sd <- 1 ##################### data <- c(-1,0,1) ##################### mean_of_data <- mean(data) assumed_sd_of_data <- sqrt(var(data)) se_data <- assumed_sd_of_data/sqrt(3) step <- 2 iter <- 10000 X<- Y <- x_prior <- y_prior<- x_likelihood <- y_likelihood <- R <-rep(0,iter) X[1] <- round(runif(1,1,5),0) for (i in 1:iter) { if (i > 1){ X[i] <- ifelse(runif(1) < R[i-1], Y[i-1], X[i-1]) } #v <- (runif(1)-0.5)*step v <- runif(1,-1,1) Y[i] <- X[i]+v x_prior[i] <- dnorm(X[i],prior_mean,prior_sd) x_likelihood[i] <- dnorm(data[1],X[i],assumed_sd_of_data)* dnorm(data[2],X[i],assumed_sd_of_data)* dnorm(data[3],X[i],assumed_sd_of_data) y_prior[i] <- dnorm(Y[i],prior_mean,prior_sd) y_likelihood[i] <- dnorm(data[1],Y[i],assumed_sd_of_data)* dnorm(data[2],Y[i],assumed_sd_of_data)* dnorm(data[3],Y[i],assumed_sd_of_data) R[i] <- exp(log(y_prior[i]) + log(y_likelihood[i])- log(x_prior[i]) - log(x_likelihood[i])) #print(c(X[i],Y[i],x_prior[i],x_likelihood[i],y_prior[i],y_likelihood[i],R[i])) } par(mfrow=c(1,2)) plot(1:iter,X,ylim=c(-4,4),type="l") abline(h=0,lty=2) mean(X[1500:3000]) sd(X[1500:3000]) plot(density(X),xlim=c(-3,5)) lines(density(rnorm(20000,prior_mean,prior_sd)),col="red") abline(v=mean(X),lty=2)