hmmheadway<-function(niter,prior,y)  #############
{n<-length(y)
 cv<-ifelse(y<mean(y),1,2)                    # Initialise component memberships.
 m<-  matrix(nrow=2,ncol=2)
 alpha<-c(1,4)                   ###################
 beta<-matrix(nrow=niter,ncol=2) ###################
 pi<-  matrix(nrow=niter,ncol=2)
 mean<-matrix(nrow=niter,ncol=2)
 proportion<-numeric(niter)
 piprop<-numeric(2)
 piold <-c(0.5,0.5)
 for (iter in 1:niter)
     {ns<-table(cv[1:(n-1)],cv[2:n])                                        # Count the transitions.
      piprop[1]<-rbeta(1,prior$a[1]+ns[2,1],prior$b[1]+ns[2,2])             # Proposal for pi_1.
      piprop[2]<-rbeta(1,prior$a[2]+ns[1,2],prior$b[2]+ns[1,1])             # Proposal for pi_2.
      Pprop<-piprop[cv[1]]/sum(piprop)                                      # Stationary probability.
      Pold <-piold[cv[1]]/sum(piold)                                        # Stationary probability.
      A<-min(1,Pprop/Pold)                                                  # Acceptance probability.
      U<-runif(1)
      if (U<A)                                                              # M-H sampling for pi.
         {pi[iter,]<-piprop
          piold<-piprop
          }
      else
         {pi[iter,]<-piold
          }
      m[1,1]<-1-pi[iter,2]                                                  # Transition matrix.
      m[1,2]<-pi[iter,1]
      m[2,1]<-pi[iter,2]
      m[2,2]<-1-pi[iter,1]
      for (comp in 1:2)                                                     # Sample other parameters.
          {nc<-sum(cv==comp)
           ycomp<-y[cv==comp]
           ysum<-sum(ycomp)                              ###################
           a<-prior$abeta[comp]+nc*alpha[comp]          ###################
           b<-prior$bbeta[comp]+ysum                    ###################
           beta[iter,comp]<-rgamma(1,a,b)                ##################
           mean[iter,comp]<-alpha[comp]/beta[iter,comp]  ##################  Mean headway.
           }
      P<-pi[iter,]/sum(pi[iter,])                                        # Stationary probabilities.
      pc<-P*m[cv[2],]                                                    # "Prior probs" for cv_1.
      for (comp in 1:2)
          {pc[comp]<-pc[comp]*dgamma(y[1],alpha[comp],beta[iter,comp])   # "Prior times likelihood". ##########
           }
      pc<-pc/sum(pc)                                                     # Normalise.
      cv[1]<-2-rbinom(1,1,pc[1])                                         # Sample cv_1.
      for (t in 2:(n-1))
          {pc<-m[,cv[t-1]]*m[cv[t+1],]                                   # "Prior probs" for cv_t.
           for (comp in 1:2)
               {pc[comp]<-pc[comp]*dgamma(y[t],alpha[comp],beta[iter,comp])  # "Prior times likelihood". #########
                }
           pc<-pc/sum(pc)                                                # Normalise.
           cv[t]<-2-rbinom(1,1,pc[1])                                    # Sample cv_t.
           }
       pc<-m[,cv[n-1]]                                                   # "Prior probs" for cv_n.
       for (comp in 1:2)
           {pc[comp]<-pc[comp]*dgamma(y[n],alpha[comp],beta[iter,comp])  # "Prior times likelihood". #########
            }
       pc<-pc/sum(pc)                                                    # Normalise.
       cv[n]<-2-rbinom(1,1,pc[1])                                        # Sample cv_n.
       proportion[iter]<-sum(cv==1)/n                                    # Proportion in component 1.
       
       }
out<-list(beta=beta,pi=pi,proportion=proportion,mean=mean)
out
 }