Gibbs sampling によるクラスター法

library(dplyr);library(mvtnorm)

n=200;m=3;dim=m-1

K=2

X=rbind(matrix(rnorm(n,mean=1,sd=1),ncol=m-1),matrix(rnorm(n,mean=5,sd=1),ncol=m-1))

mu0=rep(0,m-1)

sigma0=cov(X)

mu=rmvnorm(K,mean=mu0,sigma=sigma0)

pi0=rep(1/K,K)

times=40

for(k in 1:times){

z_mat=array(0,dim=c(n,K))

for(j in 1:nrow(X)){

x=X[j,] 

prob_vec=c()

for(i in 1:K){

prob_vec=c(prob_vec,dmvnorm(x,mean=mu[i,]))  

}

vec=prob_vec/sum(prob_vec)

z_mat[j,]=rmultinom(1,1,prob=vec)

}

n_vec=apply(z_mat,2,sum)

sign=apply(t(t(z_mat)*c(1:K)),1,sum)

X_data=data.frame(X,sign=sign)

mu_vec=array(0,dim=c(K,dim))

for(i in 1:K){

X_sub=X_data %>% filter(sign==i)  

X_sub=as.matrix(X_sub[,1:(ncol(X_sub)-1)])

ave_x=apply(X_sub,2,mean)

mu_vec[i,]=rmvnorm(1,mean=(n_vec[i]/(n_vec[i]+1))*ave_x,sigma=diag(1,dim)/(n_vec[i]+1)) 

}

mu=mu_vec

print(mu_vec)

}

library(dplyr);library(mvtnorm)

n=2000;m=7;dim=m-1

K=4

X=rbind(matrix(c(rnorm((n/2)*(m-1),mean=1,sd=1),rnorm((n/2)*(m-1),mean=5,sd=1)),ncol=m-1,nrow=n))

mu0=rep(1,m-1)

sigma0=cov(X)

tau=1;p=1;a=1;b=1

mu=rmvnorm(K,mean=mu0,sigma=sigma0)

pi0=rep(1/K,K)

times=40

for(k in 1:times){

z_mat=array(0,dim=c(n,K))

for(j in 1:nrow(X)){

x=X[j,] 

prob_vec=c()

for(i in 1:K){

prob_vec=c(prob_vec,dmvnorm(x,mean=mu[i,]))  

}

vec=prob_vec/sum(prob_vec)

z_mat[j,]=rmultinom(1,1,prob=vec)

}

n_vec=apply(z_mat,2,sum)

sign=apply(t(t(z_mat)*c(1:K)),1,sum)

X_data=data.frame(X,sign=sign)

mu_vec=array(0,dim=c(K,dim))

ave_x_mat=array(0,dim=c(K,dim))

a=n*dim/2;b=0

for(i in 1:K){

X_sub=X_data %>% filter(sign==i)  

if(nrow(X_sub)>0){

X_sub=as.matrix(X_sub[,1:(ncol(X_sub)-1)])

ave_x=apply(X_sub,2,mean)  

mu_vec[i,]=rmvnorm(1,mean=(ave_x*n_vec[i]/(n_vec[i]+p)+mu0*p/(p+n_vec[i])),sigma=diag(tau,dim)/(n_vec[i]+p)) 

ave_x_mat[i,]=ave_x

b=b+(1/2)*sum((X_sub-ave_x)^2)+(n_vec[i]*p/(2*(n_vec[i]+p)))*sum((ave_x-mu0)^2)

}

}

mu=mu_vec

tau=rgamma(1,a,b)

print(mu_vec);print(tau);print(p)

}