#mixed latent class model
#p.81
library(dplyr);library(dummies)
data=data.frame(iris)
X=data[,!(colnames(data) %in% c("Species"))]
data[,!(colnames(data) %in% c("Species"))]=t((t(X)-c(apply(X,2,mean)))/c(apply(X,2,sd)))
data[,!(colnames(data) %in% c("Species"))]=ifelse(data[,!(colnames(data) %in% c("Species"))]>0,1,0)
data$Species=as.integer(data$Species)
K=length(unique(data$Species))
Y=dummy(data$Species)
X=array(0,dim=dim(data[,!(colnames(data) %in% c("Species"))]))
for(j in 1:nrow(X)){
val=rpois(ncol(X),sample(10:20,1))
no=sample(1:ncol(X),1,,prob=val/sum(val))
X[j,no]=1
}
k=5
hyt=array(0,dim=c(k,2))
yt=array(0,dim=c(k,2))
for(j in 1:nrow(yt)){
val=rpois(2,sample(10:20,1))
hyt[j,]=val/sum(val)
no=sample(1:2,1,prob=val/sum(val))
yt[j,no]=1
}
yt=yt[,1]
hyt_vec=c()
for(j in 1:length(yt)){
hyt_vec=c(hyt_vec,ifelse(yt[j]>0,hyt[j,2],hyt[j,1]))
}
alphas=array(0.1,dim=c(ncol(X),2))
rit=array(0.1,dim=c(ncol(X),k))
Nt=rep(1,k)
L_pre=-10^10;L=-10^10+1
ite=1000;
eta=0.001
for(l in 1:ite){
#if(L-L_pre>0.001){
L_pre=L
f=array(0,dim=c(nrow(X),k))
for(j in 1:k){
for(i in 1:nrow(f)){
value1=1/(1+exp(-alphas%*%c(1,yt[j])))
value2=1-value1
fxyt=prod((c(value1)^X[i,])*(c(value2)^(1-X[i,])))
f[i,j]=fxyt
}
}
f_xh=f%*%hyt_vec
L=sum(log(f_xh))
for(t in 1:k){
for(i in 1:ncol(X)){
rit[i,t]=hyt_vec[t]*sum(X[,i]*f[,t]/f_xh)
}
}
for(t in 1:k){
Nt[t]=hyt_vec[t]*sum(f[,t]/f_xh)
}
for(j in 1:2){
for(i in 1:ncol(X)){
pis=c()
for(t in 1:k){
pis=c(pis,1/(1+exp(-sum(alphas[i,]*c(1,yt[t])))))
}
alphas[i,j]=alphas[i,j]+eta*sum(yt*(rit[i,]-Nt*pis))
}
}
print(L)
#}else{
#print("stop")
#}
}
#p.137 latent class model
library(dplyr);library(dummies)
data=data.frame(iris)
X=data[,!(colnames(data) %in% c("Species"))]
data[,!(colnames(data) %in% c("Species"))]=t((t(X)-c(apply(X,2,mean)))/c(apply(X,2,sd)))
data[,!(colnames(data) %in% c("Species"))]=ifelse(data[,!(colnames(data) %in% c("Species"))]>0,1,0)
data$Species=as.integer(data$Species)
nos=data %>% group_by(Species) %>% summarise(n=n())
nos=50
K=length(unique(data$Species))
Y=dummy(data$Species)
X=data[,!(colnames(data) %in% c("Species"))]
eta=rep(1,K);
pis=array(0.1,dim=c(ncol(X),K));
g=function(x){
return(exp(-sum(abs(x))))
}
ite=100
cost=c()
for(j in 1:ite){
mat=array(0,dim=c(nos,K))
pis_pre=pis;eta_pre=eta
for(i in 1:K){
dat=data %>% filter(Species==i) %>% select(-Species)
vec=c()
for(k in 1:nrow(dat)){
vec=c(vec,eta[i]*g(dat[k,]))
}
mat[,i]=vec
}
f=apply(mat,1,sum)
h=mat/c(f)
eta=apply(h,2,mean)
pis=array(0,dim=c(ncol(X),K))
for(i in 1:K){
dat=data %>% filter(Species==i) %>% select(-Species)
val=apply(dat*c(h[,i]),2,sum)/sum(c(h[,i]))
pis[,i]=c(val)
}
log_lik=c()
for(k in 1:nos){
likelivec=c()
for(i in 1:K){
dat=data %>% filter(Species==i) %>% select(-Species)
p=pis[,i]^(dat[k,]);q=(1-pis[,i])^(1-dat[k,])
likelihood=prod(p*q)*eta[i];
likelivec=c(likelivec,sum((as.numeric(likelihood))))
}
log_lik=c(log_lik,log(sum(likelivec)))
}
print(sum(as.numeric(log_lik)))
cost=c(cost,sum(as.numeric(log_lik)))
#print(sum(abs(eta_pre-eta))+sum(abs(pis_pre-pis)))
}
plot(cost)
#mixed latent class model p.167
#latent variable models and factor analysis(Bartholomew and Martin)
library(dplyr)
p=sample(10:50,2,replace = T)
K=10
n=100
X1=data.frame(class=0,t(rep(1,p[1])))
for(j in 1:K){
X1_sub=data.frame(class=j,matrix(rpois(n*p[1],sample(1:10,1)),ncol=p[1]))
X1=rbind(X1,X1_sub)
}
X1=tail(X1,nrow(X1)-1)
X2=data.frame(class=0,t(rep(1,p[2])))
for(j in 1:K){
X2_sub=data.frame(class=j,matrix(sample(0:1,n*p[2],replace=T),ncol=p[2]))
X2=rbind(X2,X2_sub)
}
X2=tail(X2,nrow(X2)-1)
g=function(x){
return((exp(-abs(x))))
}
eta=rep(1,K);pis=array(1,dim=c(p[2],K))
mu=array(1,dim=c(p[1],K))
sigma=rep(1,p[1])
ite=100
for(s in 1:ite){
variance1=rep(0,p[1]);variance2=0
for(j in 1:K){
X1_sub=X1 %>% filter(class==j) %>% select(-class)
X2_sub=X2 %>% filter(class==j) %>% select(-class)
val1=apply(g(X1_sub),1,sum);val2=apply(g(X2_sub),1,sum)
h=(eta[j]*(val1*val2))
eta[j]=sum(h)/n
pis_vec=apply(X2_sub*h,2,sum)/(n*eta[j])
pis[,j]=c(pis_vec)
mu_vec=apply(X1_sub*h,2,sum)/(n*eta[j])
mu[,j]=c(mu_vec)
variance1=variance1+c(apply((t(t(X1_sub)-mu_vec)^2)*h,2,sum))
variance2=variance2+sum(h)
}
sigma_pre=sigma
sigma=variance1/variance2
print(sum(abs(sigma-sigma_pre)))
}