gibbs2=function(n=1000 , m=10 , a = 1, b= 1, start1=5 , start2=0.5 )
 {
hbreaks <- c(0:(m +1)) - 0.5
xpos <- c(0:m)
ypos <- c(0:m)
for (i in 1:(m+1))
{
ypos[i] <- choose(m, xpos[i])*beta(xpos[i]+a, m-xpos[i]+b)/beta(a,b)
}
x <- seq(0,1,length=100)
y <- dbeta(x,a,b)

  x1 <- c(1:n)
  x2 <- c(1:n)
  x1[1] <- start1
  x2[1] <- start2
  par( mfrow=c( 3,2 ) )
  for ( i in 2:n )
   {
    x1[i] <- rbinom( 1 , m, x2[i-1])
    x2[i] <- rbeta(1, a + x1[i], m - x1[i] + b)
}
  plot( x1 , type="l" , xlab="Iteration" ,xlim=c(1,n),ylab="x1",main="Trace plot for x1")
  plot( x2 , type="l" , xlab="Iteration",xlim=c(1,n),ylab="x2", main="Trace plot for x2" )
  plot( x1 , x2 , type="p" , pch="." ,cex=2,xlab="x1",ylab="x2", main="Plot of each (x1,x2)")
  plot( x1 , x2 , type="l" ,col="red",xlab="x1",ylab="x2",main = "Component-by-component updating")
hist(x1,prob=T,breaks=hbreaks,xlab="x1",main="Histogram of x1 with target density")
lines(xpos,ypos)
hist(x2,prob=T,xlab="x2",main="Histogram of x2 with target density")
lines(x,y,col="black")
}