####################MCMC for GMM with one dimensions and two components########

#### libraries ########
library(stats)
library(MASS)
library(mvtnorm)
library(MCMCpack)
#####################

num = 82  ##number of obervations
numOfk=6  ##number of components
perfactorial<-factorial(numOfk);



#x<- NULL  # for real data

x<-c(
9.172, 9.350,  9.483 , 9.558 , 9.775, 10.227, 10.406, 16.084, 16.170, 18.419,
18.552, 18.600, 18.927, 19.052, 19.070, 19.330, 19.343, 19.349, 19.440, 19.473,
19.529,19.541, 19.547, 19.663, 19.846, 19.856, 19.863, 19.914, 19.918, 19.973,
19.989, 20.166, 20.175, 20.179, 20.196, 20.215,20.221, 20.415, 20.629, 20.795,
20.821, 20.846, 20.875, 20.986, 21.137, 21.492, 21.701, 21.814, 21.921, 21.960,
22.185, 22.209, 22.242, 22.249, 22.314, 22.374, 22.495, 22.746, 22.747, 22.888,
22.914, 23.206, 23.241, 23.263, 23.484, 23.538, 23.542, 23.666, 23.706, 23.711,
24.129, 24.285, 24.289, 24.366, 24.717, 24.990, 25.633, 26.960, 26.995, 32.065,
32.789, 34.279
)

#z<- NULL  # for component identity




mlist<-c(runif(1), runif(1), runif(1), runif(1), runif(1), runif(1))

##sigma = 1;

zlist = rep(0, num);


alpha <- 0.001
beta  <- 0.001

sigma  <- 100
pmean <-6
pvar<-100

mu1list<-NULL;
mu2list<-NULL;
mu3list<-NULL;
mu4list<-NULL;
mu5list<-NULL;
mu6list<-NULL;
w1list <- NULL;
w2list <- NULL;
w3list <- NULL;
w4list <- NULL;
w5list<-NULL;
w6list<-NULL;

sigmalist<-NULL;
w<-c(1/6,1/6,1/6,1/6,1/6,1/6);


M=500
B=100

for( k in 1:M){
       #update the zlist
	prob1 = (w[1])*dnorm(x, mlist[1] , sqrt(sigma)) 
        prob2 = (w[2])*dnorm(x, mlist[2] , sqrt(sigma))
        prob3 = (w[3])*dnorm(x, mlist[3] , sqrt(sigma))
        prob4 = (w[4])*dnorm(x, mlist[4] , sqrt(sigma))
        prob5 = (w[5])*dnorm(x, mlist[5] , sqrt(sigma))    
        prob4 = (w[6])*dnorm(x, mlist[6] , sqrt(sigma))         

        sumprob = prob1+prob2+prob3+prob4+prob5+prob6 
        
       #sampling and update zlist
        for (i in 1:num)   zlist[i] = sample(1:6,size=1,
                   prob=c(prob1[i]/sumprob[i],prob2[i]/sumprob[i], prob3[i]/sumprob[i], prob4[i]/sumprob[i],prob5[i]/sumprob[i], prob6[i]/sumprob[i]))
 
	#update the sigma
          sigma = 1/rgamma(1,alpha+num/2,beta+sum((x-mlist[zlist])^2)/2)
          sigmalist<-c(sigmalist, sigma);
	
 	#update mu1

	numOfgroup1 = sum( zlist==1 );
        numOfgroup2 = sum( zlist==2 );
        numOfgroup3 = sum( zlist==3 );
        numOfgroup4 = sum( zlist==4 );
        numOfgroup5 = sum( zlist==5 );
        numOfgroup6 = sum( zlist==6 );
        
        var  = 1/(1/pvar+ numOfgroup1/sigma)
        mean = var*( pmean/pvar+sum(x[zlist==1])/sigma)
        
        mlist[1] = rnorm(1,mean,sqrt(var))


        var  = 1/(1/pvar+ numOfgroup2/sigma)
        mean = var*( pmean/pvar+sum(x[zlist==2])/sigma)
        
        mlist[2] = rnorm(1,mean,sqrt(var))
	
	var  = 1/(1/pvar+ numOfgroup3/sigma)
        mean = var*( pmean/pvar+sum(x[zlist==3])/sigma)
        
        mlist[3] = rnorm(1,mean,sqrt(var))
 
        var  = 1/(1/pvar+ numOfgroup4/sigma)
        mean = var*( pmean/pvar+sum(x[zlist==4])/sigma)
        
        mlist[4] = rnorm(1,mean,sqrt(var)) 

	var  = 1/(1/pvar+ numOfgroup5/sigma)
        mean = var*( pmean/pvar+sum(x[zlist==5])/sigma)
        
        mlist[5] = rnorm(1,mean,sqrt(var)) 

        var  = 1/(1/pvar+ numOfgroup6/sigma)
        mean = var*( pmean/pvar+sum(x[zlist==6])/sigma)
        
        mlist[6] = rnorm(1,mean,sqrt(var)) 

        #mu2list <- c(mu2list, mlist[2]);
        

	###Random Permutation    
	 swith <- sample(1:6);
	 temp1 <- mlist[swith[1]];
         temp2 <- mlist[swith[2]];
         temp3 <- mlist[swith[3]];
         temp4 <- mlist[swith[4]];
         temp5 <- mlist[swith[5]];
         temp6 <- mlist[swith[6]];
        
         mlist[1]<-temp1;
         mlist[2]<-temp2;
         mlist[3]<-temp3;
         mlist[4]<-temp4;						
         mlist[5]<-temp5;
         mlist[6]<-temp6;  


				
        mu1list<-c(mu1list, mlist[1])
      	mu2list <- c(mu2list, mlist[2]);
        mu3list<-c(mu3list, mlist[3])
      	mu4list <- c(mu4list, mlist[4]);
        mu5list<-c(mu5list, mlist[5])
      	mu6list <- c(mu6list, mlist[6]);
	

        #update w
	
        w = rdirichlet( 1, c(1+numOfgroup1, 1+numOfgroup2, 1+numOfgroup3, 1+numOfgroup4, 1+numOfgroup5, 1+numOfgroup6) )
    	w1list<-c(w1list, w[1]);
        w2list<-c(w2list, w[2]);
	w3list<-c(w3list, w[3]);
	w4list<-c(w4list, w[4]);
	w5list<-c(w5list, w[5]);
	w6list<-c(w6list, w[6]);

     

}

cat("mean of mu1 is: " , mean(mu1list[(B+1):M]) , "\n");
cat("mean of mu2 is: " , mean(mu2list[(B+1):M]) , "\n");
cat("mean of mu3 is: " , mean(mu3list[(B+1):M]) , "\n");
cat("mean of mu4 is: " , mean(mu4list[(B+1):M]) , "\n");
cat("mean of mu5 is: " , mean(mu5list[(B+1):M]) , "\n");
cat("mean of mu6 is: " , mean(mu6list[(B+1):M]) , "\n");