##### Below is the Splus code for the last example in Ch5, part I.



K_5000   ### K is the number of samples I want to draw from the posterior distribution

x_c(0:1000)/1000  	### x represents theta. The interval is divided into 1000 pieces, each of length 1/1000.
			### Here I am let x[1]=0, x[2]=1/1000, ...., x[1000]=999/1000, x[1001]=1.
dx_rep(0,1001)        	### dx represents the density at x. Here I initialize dx as a zero vector of length 1001. 
			 

dx[1]_0			
dx[1001]_0		### Here I calculate dx at the end point x=0 and x=1.

dx[c(2:500)]_exp(61*log(x[(2:500)])+40*log(1-x[c(2:500)])) ## calculate dx for x<=0.5. Here log = ln.

dx[c(501:1000)]_exp(60*log(x[c(501:1000)])+41*log(1-x[c(501:1000)]))  ## calculate dx for x>0.5

dx_dx/sum(dx)*1000      ### This is the normalization part in order to make a discrete probability distribution.
			### dx now is a real discrete probability function.

dy_rep(0,1001)          ### dy will represent the cumulative distribution function (cdf) for the discrete approximation.
			### This is initialization of dy.

for (i in 1:1001)	### Calculate dy. dy[i] represents the cdf at x[i]. or dy[i] = (dx[1]+dx[2]+...dx[i])/1000.	
{
  dy[i]_sum(dx[c(1:i)])/1000
}

z_runif(K,0,1)		### generate K uniform samples from [0,1].

ind_rep(1,K)            ### for each uniform sample, ind will represent the corresponding (where the uniform 
			### sample falls) in your class note
for (i in 1:1000)
{
  ind[z>dy[i]]_i
}  

hist(x[ind],breaks=c(0:90)/90, probability=T, xlab="") ### x[ind] = x[j] is the corresponding sample from the discrete approximation.
 							### hist() is to draw the histgram