#混合指数分布によるクラスタリング
#mixed exponential Intervals p.29
n=1000
k=4
X=array(0,dim=c(k,n))
param=c()
for(j in 1:k){
val=sample(5:50,1)
param=c(param,val)
X[j,]=abs(rpois(n,val))
}
X=as.numeric(X)
lam=1/sample(5:10,k);p=rep(1/k,k)
ite=1000
for(l in 1:ite){
q=array(0,dim=c(k,length(X)))
for(i in 1:k){
for(j in 1:length(X)){
q[i,j]=p[i]*dexp(X[j],1/lam[i])
}
}
q=t(t(q)/c(apply(q,2,sum)))
p=apply(q,1,mean)
lam=1/(q%*%X/c(apply(q,1,sum)))
q[q==0]=10^(-20)
#print(sum(log(q)))
print(lam)
}
#混合正規分布によるクラスタリング
library(dplyr)
N=10000;
class=4
pi=rep(0.25,class)
mu=c(-5,5,8,3);sig=c(1,5,3,4)
#真の出身（正解データ）
attr=sample(1:class,N,replace = T,prob=pi)
x=rep(-99,N)
for(j in 1:class){
x[which(attr==j)]<-rnorm(length(which(attr==j)),mu[j],sig[j])
}
#initial conditions
mu=sample(c(-20:20),class);sig=sample(c(1:50),class)
pi=c(0.2,0.6,0.1,0.1)
times=100
for(j in 1:times){
s_piN=array(0,dim=c(class,N))
for(i in 1:class){
s_piN[i,]=pi[i]*dnorm(x,mu[i],sig[i])
}
qn=c()
for(i in 1:class){
qn=c(qn,sum(s_piN[i,]/apply(s_piN,2,sum))/N)
}
pi=qn
mu=c()
sig=c()
for(i in 1:class){
mu_val=sum(x*s_piN[i,]/apply(s_piN,2,sum))/(N*pi[i])
mu=c(mu,mu_val)
sig=c(sig,sqrt(sum((s_piN[i,]/apply(s_piN,2,sum))*((x-mu_val)*(x-mu_val)))/(N*pi[i])))
}
}
N=1000
pi0=0.5;pi1=0.5
mu0=5;mu1=-5
sig0=1;sig1=5
attr=sample(0:1,N,replace=T,prob=c(pi0,pi1))
x=rep(-99,N)
x[which(attr==0)]=rnorm(length(which(attr==0)),mu0,sig0)
x[which(attr==1)]=rnorm(length(which(attr==1)),mu1,sig1)
#initial conditions
mu0=1;mu1=1;sig0=10;sig1=20;pi0=0.1;pi1=0.9
iteration=20
for(i in 1:iteration){
piN0=pi0*dnorm(x,mu0,sig0);piN1=pi1*dnorm(x,mu1,sig1)
qn0=piN0/(piN0+piN1);qn1=piN1/(piN0+piN1)
pi0=sum(qn0)/N;pi1=sum(qn1)/N
mu0=sum(qn0*x)/(N*pi0);mu1=sum(qn1*x)/(N*pi1)
sig0=sqrt(sum(qn0*(x-mu0)*(x-mu0))/(N*pi0))
sig1=sqrt(sum(qn1*(x-mu1)*(x-mu1))/(N*pi1))
}