RBM

```{r}

#制限ボルツマンマシン

#CD method 未完

library(dplyr)

data=data.frame(iris)

data$Species=as.integer(data$Species)

X=data[,colnames(data)!=c("Species")]

K=50;N=nrow(X)

a=sample(1:10,ncol(X),replace=T)

w=array(2,dim=c(length(a),K))

V_mat=t((t(X)-apply(X,2,mean))/apply(X,2,sd))

b=sample(10:20,K,replace = T);v=sample(0:1,length(a),replace = T)

h=sample(0:1,length(b),replace=T)


ite=10000

eta=0.01

v=V_mat[sample(1:nrow(V_mat),1),];h_vec=h

v_vec=V_mat[sample(1:nrow(V_mat),1),]

p_h_0=exp(b+apply(w*c(v_vec),2,sum))/(1+exp(b+apply(w*c(v_vec),2,sum)))

p_v_0=exp(a+apply(t(w)*c(h_vec),2,sum))/(1+exp(a+apply(t(w)*c(h_vec),2,sum)))


cost=c();

for(l in 1:ite){
  
pre_a=a;pre_b=b;pre_w=w  
  
p_h=exp(b+apply(w*c(v_vec),2,sum))/(1+exp(b+apply(w*c(v_vec),2,sum)))  
  
p_v=exp(a+apply(t(w)*c(h_vec),2,sum))/(1+exp(a+apply(t(w)*c(h_vec),2,sum)))

#v_vec=ifelse(sample(seq(0,1,0.001),length(a),replace=T)<p_v,1,0)

#h_vec=ifelse(sample(seq(0,1,0.001),length(b),replace=T)<p_h,1,0)

v_vec=c()

for(j in 1:length(a)){
  
v_vec=c(v_vec,sample(0:1,1,prob=c(1-p_v[j],p_v[j])))  
  
}

h_vec=c()

for(j in 1:length(b)){
  
h_vec=c(h_vec,sample(0:1,1,prob=c(1-p_h[j],p_h[j])))  
  
}

for(i in 1:length(a)){
for(j in 1:length(b)){  

w[i,j]=w[i,j]+eta*(v[i]*p_h_0[j]-v_vec[i]*p_h[j])    
  
}
}

for(i in 1:length(a)){
  
a[i]=a[i]+eta*(v[i]-v_vec[i])

}

for(j in 1:length(b)){
  
b[j]=b[j]+eta*(p_h_0[j]-p_h[j])  
  
}

#print(sum(abs(pre_w-w))+sum(abs(pre_a-a))+sum(abs(pre_b-b)))

cost=c(cost,sum(abs(t(t(v))%*%t(p_h_0)-t(t(v_vec))%*%t(p_h)))+sum(abs(v-v_vec))+sum(abs(p_h-p_h_0)))

print(sum(abs(t(t(v))%*%t(p_h_0)-t(t(v_vec))%*%t(p_h)))+sum(abs(v-v_vec))+sum(abs(p_h-p_h_0)))

}  
  

plot(cost)  



```

```{r}

#可視変数vが連続

library(dplyr)

data=data.frame(iris)

data$Species=as.integer(data$Species)

X=data[,colnames(data)!=c("Species")]

K=500;N=nrow(X)

a=sample(1:10,ncol(X),replace=T)

w=array(2,dim=c(length(a),K))

#V_mat=t((t(X)-apply(X,2,mean))/apply(X,2,sd))

V_mat=as.matrix(X)/100

b=sample(10:20,K,replace = T);v=sample(0:1,length(a),replace = T)

h=sample(0:1,length(b),replace=T)


ite=10000

eta=0.1;M=1000

cost=c();

for(l in 1:ite){
  
pre_a=a;pre_b=b;pre_w=w  

v_vec1=array(0,dim=c(M,length(a)))

for(i in 1:M){
  
v_vec1[i,]=as.numeric(V_mat[sample(1:nrow(V_mat),1),])

}

V_E_p_h=array(0,dim=c(nrow(V_mat),K))

for(i in 1:nrow(V_mat)){

v_vec=V_mat[i,]  
    
p_h=exp(b+apply(w*c(v_vec),2,sum))/(1+exp(b+apply(w*c(v_vec),2,sum)))  

V_E_p_h[i,]=p_h
  
}


E_p_h=array(0,dim=c(M,K));E_p_hv=c()

for(i in 1:M){
  
p_h=exp(b+apply(w*c(v_vec1[i,]),2,sum))/(1+exp(b+apply(w*c(v_vec1[i,]),2,sum)))  
 
E_p_h[i,]=p_h

}

cost_vec=array(0,dim=c(length(a),length(b)))

for(i in 1:length(a)){
for(j in 1:length(b)){  

w[i,j]=w[i,j]+eta*(sum(V_mat[,i]*V_E_p_h[,j])-sum(v_vec1[,i]*E_p_h[,j]))

cost_vec[i,j]=(sum(V_mat[,i]*V_E_p_h[,j])-sum(v_vec1[,i]*E_p_h[,j]))

}
}

b=b+eta*(apply(V_E_p_h,2,mean)-apply(E_p_h,2,mean))

a=a+eta*(apply(V_mat,2,mean)-apply(v_vec1,2,mean))

print(sum(abs(apply(V_E_p_h,2,mean)-apply(E_p_h,2,mean)))+sum(abs(apply(V_mat,2,mean)-apply(v_vec1,2,mean)))+sum(abs(cost_vec)))

}  
  




```

```{r}



#可視変数vが連続 改変

library(dplyr)

data=data.frame(iris)

data$Species=as.integer(data$Species)

X=data[,colnames(data)!=c("Species")]

K=500;N=nrow(X)

a=sample(1:10,ncol(X),replace=T)

w=array(2,dim=c(length(a),K))

#V_mat=t((t(X)-apply(X,2,mean))/apply(X,2,sd))

V_mat=as.matrix(X)/100

b=sample(10:20,K,replace = T);v=sample(0:1,length(a),replace = T)

h=sample(0:1,length(b),replace=T)


ite=100

eta=0.1;M=1000

cost=c();

for(l in 1:ite){

pre_a=a;pre_b=b;pre_w=w  

#v_vec1=array(0,dim=c(M,length(a)))

v_vec1=V_mat[sample(1:nrow(V_mat),M,replace=T),]

#for(i in 1:M){

#v_vec1[i,]=as.numeric(V_mat[sample(1:nrow(V_mat),1),])

#}

#V_E_p_h=array(0,dim=c(nrow(V_mat),K))

V_E_p_h=exp(t(t(w)%*%t(V_mat))+c(b))/(1+exp(t(t(w)%*%t(V_mat))+c(b)))

#for(i in 1:nrow(V_mat)){

#v_vec=V_mat[i,]  

#p_h=exp(b+apply(w*c(v_vec),2,sum))/(1+exp(b+apply(w*c(v_vec),2,sum)))  

#V_E_p_h[i,]=p_h

#}


#E_p_h=array(0,dim=c(M,K));E_p_hv=c()

E_p_h=exp(t(w)%*%t(v_vec1)+c(b))/(1+exp(t(w)%*%t(v_vec1)+c(b)))

#for(i in 1:M){

#p_h=exp(b+apply(w*c(v_vec1[i,]),2,sum))/(1+exp(b+apply(w*c(v_vec1[i,]),2,sum)))  
#E_p_h[i,]=p_h

#}

cost_vec=t(V_mat)%*%V_E_p_h-t(v_vec1)%*%t(E_p_h)

w=w+eta*cost_vec

b=b+eta*(apply(V_E_p_h,2,mean)-apply(E_p_h,2,mean))

a=a+eta*(apply(V_mat,2,mean)-apply(v_vec1,2,mean))

print(sum(abs(apply(V_E_p_h,2,mean)-apply(E_p_h,2,mean)))+sum(abs(apply(V_mat,2,mean)-apply(v_vec1,2,mean)))+sum(abs(cost_vec)))

}  





  
```