####################MCMC for GMM with one dimensions and two components######## #### libraries ######## library(stats) library(MASS) library(mvtnorm) library(MCMCpack) ##################### num = 82 ##number of obervations numOfk=3 ##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 ########The MCMC start from here ############## mlist<-c(runif(1), runif(1), runif(1)) ##sigma = 1; zlist = rep(0, num); alpha <- 0.001 beta <- 0.001 sigma <- c(100, 100, 100); pmean <-6 pvar<-100 mu1list<-NULL; mu2list<-NULL; mu3list<-NULL; mu4list<-NULL; w1list <- NULL; w2list <- NULL; w3list <- NULL; w4list <- NULL; sigmalist1<-NULL; sigmalist2<-NULL; sigmalist3<-NULL; w<-c(1/3,1/3,1/3); M=10000 B=1000 for( k in 1:M){ #update the zlist prob1 = (w[1])*dnorm(x, mlist[1] , sqrt(sigma[1])) prob2 = (w[2])*dnorm(x, mlist[2] , sqrt(sigma[2])) prob3 = (w[3])*dnorm(x, mlist[3] , sqrt(sigma[3])) # prob4 = (w[4])*dnorm(x, mlist[4] , sqrt(sigma)) sumprob = prob1+prob2+prob3 #sampling and update zlist for (i in 1:num) zlist[i] = sample(1:numOfk,size=1, prob=c(prob1[i]/sumprob[i],prob2[i]/sumprob[i], prob3[i]/sumprob[i])) numOfgroup1 = sum( zlist==1 ); numOfgroup2 = sum( zlist==2 ); numOfgroup3 = sum( zlist==3 ); #numOfgroup4 = sum( zlist==4 ); #update the sigma sigma[1] = 1/rgamma(1,(1+alpha+numOfgroup1)/2,beta+sum((x[zlist==1]-mlist[1])^2)/2) sigma[2] = 1/rgamma(1,(1+alpha+numOfgroup2)/2,beta+sum((x[zlist==2]-mlist[2])^2)/2) sigma[3] = 1/rgamma(1,(1+alpha+numOfgroup3)/2,beta+sum((x[zlist==3]-mlist[3])^2)/2) sigmalist1<-c(sigmalist1, sigma[1]); sigmalist2<-c(sigmalist2, sigma[2]); sigmalist3<-c(sigmalist3, sigma[3]); #update mu1 var = 1/(1/pvar+ numOfgroup1/sigma[1]) mean = var*( pmean/pvar+sum(x[zlist==1])/sigma[1]) mlist[1] = rnorm(1,mean,sqrt(var)) var = 1/(1/pvar+ numOfgroup2/sigma[2]) mean = var*( pmean/pvar+sum(x[zlist==2])/sigma[2]) mlist[2] = rnorm(1,mean,sqrt(var)) var = 1/(1/pvar+ numOfgroup3/sigma[3]) mean = var*( pmean/pvar+sum(x[zlist==3])/sigma[3]) mlist[3] = rnorm(1,mean,sqrt(var)) ###Random Permutation swith <- sample(1:4); temp1 <- mlist[swith[1]]; temp2 <- mlist[swith[2]]; temp3 <- mlist[swith[3]]; mlist[1]<-temp1; mlist[2]<-temp2; mlist[3]<-temp3; mu1list<-c(mu1list, mlist[1]) mu2list <- c(mu2list, mlist[2]); mu3list<-c(mu3list, mlist[3]) #mu4list <- c(mu4list, mlist[4]); #update w w = rdirichlet( 1, c(1+numOfgroup1, 1+numOfgroup2, 1+numOfgroup3) ) w1list<-c(w1list, w[1]); w2list<-c(w2list, w[2]); w3list<-c(w3list, w[3]); } ###Q matrix #Pij<-matrix(0, num, 4); pij<-array(1/3, c(M, 6, num,3)); meanlist<-cbind(mu1list, mu2list, mu3list); sigmalist<-cbind(sigmalist1, sigmalist2, sigmalist3); pilist <- cbind(w1list, w2list, w3list); #####all the perm<-rbind(c(1,2,3), c(1,3,2), c(2,1,3),c(2,3,1),c(3,1,2),c(3,2,1)); pmatrix=function(i,j,t,m){ a=0 for (k in 1:3) { a = a+pilist[m,k] * dnorm(x[i], meanlist[m,perm[t,k]] , sqrt(sigmalist[m,perm[t,k]]) ); } pilist[m,j]*dnorm(x[i], meanlist[m,perm[t,j]],sqrt(sigmalist[m,perm[t,j]]) )/a } minpij<-function(s, t, Qij){ # sum( pij[s, t,,]*log(pij[s,t,,]/Qij)) sum( exp(pij[s, t,,]) + pij[s,t,,]/Qij ); } ######compute the pij matrix############################################ for (m in 1:M) { for (t in 1:6 ){ for (i in 1:num){ for (j in 1:3){ pij[m,t,i,j]<-pmatrix(i,j,t,m); } } } #end of m for loop } #end of p for loop oldperm<-rep(1,M); newperm<-rep(1,M); exi = 0; while (exi != 1 ){ ##update Q Q=matrix(0, num,3); for( m in 1:M){ Q<-pij[m,oldperm[m],,]+Q; } Q<-Q/M; ###choose permutation for( s in 1:M){ ##loop through MCMC samples minlist<-NULL; for(t in 1:6){ #for each permutation, find the pij value minlist<-c(minlist, minpij( s, t, Q) ); } mindex <- which.min(minlist); newperm[s]<-mindex; } if( sum(oldperm==newperm)>M-10 ) {exi = 1} else { oldperm = newperm; } } newmu1list<-NULL; newmu2list<-NULL; newmu3list<-NULL; newsigmalist1<-NULL; newsigmalist2<-NULL; newsigmalist3<-NULL; for( s in 1:M){ newmu1list <- c(newmu1list, meanlist[ s, perm[newperm[s],1] ] ); newmu2list <- c(newmu2list, meanlist[ s, perm[newperm[s],2] ] ); newmu3list <- c(newmu3list, meanlist[ s, perm[newperm[s],3] ] ); newsigmalist1<- c(newsigmalist1, meanlist[ s, perm[newperm[s],1] ] ); newsigmalist2 <- c(newsigmalist2, meanlist[ s, perm[newperm[s],2] ] ); newsigmalist3<- c(newsigmalist3, meanlist[ s, perm[newperm[s],3] ] ); } cat("mean of newmu1 is: " , mean(newmu1list[(B+1):M]) , "\n"); cat("mean of newmu2 is: " , mean(newmu2list[(B+1):M]) , "\n"); cat("mean of newmu3 is: " , mean(newmu3list[(B+1):M]) , "\n"); 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");