# 制限ボルツマンマシン
# 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)
# 可視変数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)))
}
# 可視変数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)))
}