set.seed(1); # save MC seed for replicablity n<- 100; y<- rnorm(n); z<- rnorm(n); x<- z + rnorm(n); my.density<- function(beta) { exp(-4*sum( (ifelse(y < x*beta, 1, 0 ) -.5)*z)^2/var(z)/n ) } # unnormalized posterior density without integration constant grids<- seq(-1,1, .0091) v<- grids eval<- apply(matrix(grids, length(grids), 1), 1, my.density) plot( grids, eval, type="l", xlab="theta", ylab="density = exp(-L(theta)" ) number.draws<- 1000; MC.chain<- function(number.draws) { start<- 0 chain<- rep(start, number.draws) for( i in 2: number.draws) { xi<- rnorm(1, mean=chain[i-1] , sd=1) # draw rej <- (my.density(xi)/my.density(chain[i-1]) ) if ( runif(1) < rej ) chain[i]<- xi else chain[i]<- chain[i-1]; } return(chain); } chain<- MC.chain(10000); plot(chain, main="Markov Chain Sequence", type="l", col=2) q.50<- median(chain, .50); # posterior median q.05<- quantile(chain, .05); # posterior .05 quantile q.95<- quantile(chain, .95); # posterior .95 quantile mm<- mean(chain); # posterior of the chain plot(density(chain), type="l", xlab="theta", ylab="Q-posterior", main="Q-Posterior for Theta", col=1) # command "density (x) " constructs a hystogram of the sample "x"; type help(density) to find out points( q.50, 0, cex=2, pch=19) points( q.95, 0, cex=2, pch=20) points( q.05, 0, cex=2, pch=20)