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)