library(diptest)
library(mixAK) 

# y is the data vector
# keep is the number of thinned MCMC iterates kept after burnin
modes=function(y,keep){
 a=min(y); b=max(y); ra=b-a; a=a-ra/4; b=b+ra/4
 gp=100; x=seq(a,b,(b-a)/(gp-1)); d=x
 d=rep(0,gp); modes=rep(0,keep)
 r=NMixMCMC(y,scale=list(shift=0,scale=1),
 prior=list(priorK="tpoisson",Kmax=20,lambda=1,xi=0.5*(a+b),ce=1/ra^2),
 nMCMC=c(burn=200,keep=keep,thin=10,info=100),keep.chains=TRUE)
 i=0
 for(j in 1:keep){
  kt=r$K[j]; w=r$w[(i+1):(i+kt)]; mu=r$mu[(i+1):(i+kt)]; si=sqrt(r$Sigma[(i+1):(i+kt)])
  for(k in 1:gp){d[k]=sum(w*dnorm(x[k],mu,si))}
  for(k in 3:gp){if(d[k-2]<d[k-1]){if(d[k-1]>d[k]){modes[j]=modes[j]+1}}}
  i=i+kt; if(modes[j]==0){modes[j]=1}
 }
 ktot=max(r$K); mc=matrix(0,ncol=2,nrow=ktot)
 for(k in 1:ktot){mc[k,1]=sum(r$K==k)/keep; mc[k,2]=sum(modes==k)/keep}
 rownames(mc)=1:max(r$K); colnames(mc)=c("P(comps)","P(modes)")
 cat("BF for >1 mode=",(19/(sum(modes==1)/(keep-sum(modes==1)))),"\n")
 print("Posterior components/modes..."); mc
}
segments=c(23,30,54,28,31,29,34,35,30,
           27,21,43,51,35,51,49,35,24,
           26,29,21,29,37,27,28,33,33,
           23,37,27,40,48,41,20,30,57)

modes(segments,2000)

hist(segments,freq=F) 
lines(density(segments),col="blue") # smoothed histogram!

dip.test(segments) # test for more than one mode
dip.test(segments,simulate.p.value=T)