metropolis.update =function(start=1, n=10, sig.q=1, mu.p=0, sig.p=1) {

  #set up vector theta[] for the samples
  theta=rep(NA,n)
  theta[1]=start

  #count the number of accepted proposals
  accept=0

  #Set up plots
  layout(matrix(c(0,1,1,0,0,1,1,0,2,2,2,2,2,2,2,2),4,4,byrow=T))
  #Selection of points
  plot (c(-5,5),c(0,0.5),type="n",main="",xlab="theta",ylab="density")
  #Target density
  curve(dnorm(x,mu.p,sqrt(sig.p)),add=T,lty=2,col="blue")
  #Proposal density
  curve(dnorm(x,theta[1],sqrt(sig.q)),add=T,col="red")
  #current proposal
  lines(c(theta[1],theta[1]),c(0,0.5),col="red")

  #Time series plot so far:
  plot(c(0,n),c(min(start,-3),max(start,3)),
       xlab="iterations",ylab="theta",type="n")
  lines(c(0,n),c(1.96,1.96),lty=2,col="blue")
  lines(c(0,n),-c(1.96,1.96),lty=2,col="blue")

  readLines(n=1)
  
  #Main loop
  for (t in 2:n) {

    #pick a candidate value theta* from N(theta[t-1], sig.q)
    theta.star = rnorm(1, theta[t-1], sqrt(sig.q))
    
    #generate a random number between 0 and 1:
    u = runif(1)

    # calculate acceptance probability:
    r = (dnorm(theta.star,mu.p,sig.p)) / (dnorm(theta[t-1],mu.p,sig.p) )
    a=min(r,1)
    

    #if u<a we accept the proposal; otherwise we reject
    if (u<a) {   
        #accept
            accept=accept+1
            theta[t] = theta.star
            cat("\ttheta*=",theta.star,"\n")
            cat("ACCEPT: theta[t]=",theta[t],"\n")
    }
    else {
        #reject
        theta[t] = theta[t-1]
        cat("\ttheta*=",theta.star,"\n")
        cat("REJECT: theta[t]=",theta[t],"\n")
        }

    # update plots
    layout(matrix(c(0,1,1,0,0,1,1,0,2,2,2,2,2,2,2,2),4,4,byrow=T))
    #Selection of points
    plot (c(-5,5),c(0,0.5),type="n",main="",xlab="theta",ylab="density")
    #Target density
    curve(dnorm(x,mu.p,sqrt(sig.p)),add=T,lty=2,col="blue")
    #Proposal density
    curve(dnorm(x,theta[t-1],sqrt(sig.q)),add=T,col="red")
    #current proposal
    lines(c(theta.star,theta.star),c(0,0.5),col="red",lty=2)
    readLines(n=1)
    #theta[t]
    lines(c(theta[t],theta[t]),c(0,0.5),col="red") 

    #Plot the density so far
    if (t>2) lines(density(theta[1:t]),col="green",lty=2)
    #Time series plot so far:
    plot(c(0,n),c(min(start,-3),max(start,3)),
         xlab="iterations",ylab="theta",type="n")
    lines(0:t,theta[1:(t+1)],col="red")
    lines(c(0,n),c(1.96,1.96),lty=2,col="blue")
    lines(c(0,n),-c(1.96,1.96),lty=2,col="blue")

       
    #Wait for keypress
    readLines(n=1)
    
  #end loop
  }


}