Latent variable Models and Factor Analysis,  David J.Bartholomew  Martin Knott

```{r}

#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")  
  
#}

}
  
```



```{r}

#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)

```

```{r}

#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)))

}


```