hmmnorm<-function(niter,prior,y)
{n<-length(y)
 cv<-ifelse(y<mean(y),1,2)                    # Initialise component memberships.
 m<-  matrix(nrow=2,ncol=2)
 mu<- matrix(nrow=niter,ncol=2)
 tau<-matrix(nrow=niter,ncol=2)
 pi<- 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)
           if (nc>0)
              {ycomp<-y[cv==comp]
               ybar<-mean(ycomp)
               s2n<-(sum(ycomp*ycomp)-nc*ybar*ybar)/nc
               ycd<-ycomp-prior$m[comp]
               r2<-sum(ycd*ycd)/nc
               k1<-prior$c[comp]+nc
               d1<-prior$d[comp]+nc
               m1<-(prior$c[comp]*prior$m[comp]+nc*ybar)/k1
               vd<-(prior$c[comp]*r2+nc*s2n)/k1
               v1<-(prior$d[comp]*prior$v[comp]+nc*vd)/d1
               tau[iter,comp]<-rgamma(1,(d1/2),(d1*v1/2))
               sd<-sqrt(1/(k1*tau[iter,comp]))
               mu[iter,comp]<-rnorm(1,m1,sd)
               }
            else
              {tau[iter,comp]<-rgamma(1,(prior$d[comp]/2),(prior$d[comp]*prior$v[comp]/2))
               sd<-sqrt(1/(prior$c[comp]*tau[iter,comp]))
               mu[iter,comp]<-rnorm(1,prior$m[comp],sd)
               }
           }
      P<-pi[iter,]/sum(pi[iter,])                                        # Stationary probabilities.
      pc<-P*m[cv[2],]                                                    # "Prior probs" for cv_1.
      sd<-numeric(2)
      for (comp in 1:2)
          {sd[comp]<-sqrt(1/tau[iter,comp])
           pc[comp]<-pc[comp]*dnorm(y[1],mu[iter,comp],sd[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]*dnorm(y[t],mu[iter,comp],sd[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]*dnorm(y[n],mu[iter,comp],sd[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(mu=mu,tau=tau,pi=pi,proportion=proportion)
out
 }