【总结】关于玻尔兹曼机(BM)、受限玻尔兹曼机(RBM)、深度玻尔兹曼机(DBM)、深度置信网络(DBN)理论总结和代码实践前言玻尔兹曼机(BM)受限玻尔兹曼机(RBM)深度玻尔兹曼机(DBM)深度置信网络(DBN)

59 阅读 0 评论 39 点赞

我是靠谱客的博主英勇百合，最近开发中收集的这篇文章主要介绍【总结】关于玻尔兹曼机(BM)、受限玻尔兹曼机(RBM)、深度玻尔兹曼机(DBM)、深度置信网络(DBN)理论总结和代码实践前言玻尔兹曼机(BM)受限玻尔兹曼机(RBM)深度玻尔兹曼机(DBM)深度置信网络(DBN)，觉得挺不错的，现在分享给大家，希望可以做个参考。

概述

近期学习总结

前言
玻尔兹曼机(BM)
- 波尔兹曼分布推导过程
- 吉布斯采样
受限玻尔兹曼机(RBM)
- 能量函数
- CD学习算法
- 代码实现受限玻尔兹曼机
深度玻尔兹曼机(DBM)
- 代码实现深度玻尔兹曼机
深度置信网络(DBN)
- 代码实现深度置信网络

前言

前段时间，团队被老师安排了一个涉及玻尔兹曼机及相关变体的任务。具体内容就是学习相关理论知识及代码实现。
所以就想着写一篇博客来总结一下繁杂的知识点（尤其是背后的数学推导公式），如果你是玻尔兹曼机及相关变体的初学者，并对这些知识非常感兴趣，那么恭喜你，发现了这篇文章。我会用非常白话、通俗易懂的方式向读者解释这些知识逻辑。
在此特别感谢王志强、董倩等研究生，对此篇博客做出的贡献。

玻尔兹曼机(BM)

玻尔兹曼机(Boltzman Machine)是一个随机动力系统，每个变量的状态都以一定的概率受到其他变量的影响。
玻尔兹曼机可以用概率无向图模型来描述，一个具有K个节点的玻尔兹曼机满足以下三个性质：

二值化。每个节点的状态值只有0和1。
一个玻尔兹曼机包括两类节点，一类是可观察的节点有N个，一类是不可观察的节点，即隐藏节点，有(K-N)个。
节点之间是全连接的。每个节点都和其他节点连接。
每两个变量之间的互相影响是对称的。这里的对称和上面无向其实是一个概念，说白了就是已知A点的状态值，那么求B的状态值和已知B的状态值，求A的状态值的影响是相等的。如果你还是没有理解这句话，碰巧你又了解过一点概率论的知识，那么你可以将上述理解为 P(B|A) = P(A|B)。

在这里插入图片描述
上图就是一个有六个节点的玻尔兹曼机。其中有三个可观察的节点，我已经标了黄色，还有三个不可观测的节点，即隐藏节点，我已经标了灰色。

波尔兹曼分布推导过程

玻尔兹曼机中，随机向量X的联合概率，也就是节点的状态值，是满足玻尔兹曼分布的。

玻尔兹曼分布是描述粒子处于特定状态下的概率，是关于状态能量E(x)与系统温度T的函数。一个粒子处于状态α的概率P(α)是关于状态能量E(x)与系统温度T的函数。

别急，我相信你读完这个定义一定是懵逼，心里已经开始呐喊，TMD啥叫特定状态啊？啥叫特定状态下的概率啊？啥是状态能量啊？…
不急，我这就用 “人话” 翻译一遍。
特定状态就是说，节点的状态值为1，还是0。x=α
特定状态就是说，节点状态值为1或者0时的概率。P(x=α)
状态能量就是说，粒子本身具有的能量。
玻尔兹曼分布就是计算P(x=α)时具体的概率函数与系统的状态能量E(x)和系统温度T的函数有关。具体表达式为：
在这里插入图片描述
通常玻尔兹曼分布还有另一个表达式，实则是上式的等价处理，如下：

E(x)为能量函数，T为系统温度，Z为配分函数，实则就是一个归一化因子。
从玻尔兹曼分布的定义，我们可以发现系统的两个不同状态的概率之比仅与系统能量有关：
$P (a) / P (b) = e x p ((E 1 - E 2) / K T)$
具体的推导公式如下图：
在这里插入图片描述

吉布斯采样

我采用的是Python下的tensooflow语言来实现的玻尔兹曼机。请注意，我们使用的版本是tensorflow2。
在讲代码之前，有必要讲一下玻尔兹曼机的训练过程。
在玻尔兹曼机中，配分函数Z通常很难计算，因此，联合概率分布P(x)一般通过马尔科夫链蒙特卡洛方法(MCMC方法)来做近似计算。
玻尔兹曼机采用了基于吉布斯采样的样本生成方法来训练的。

吉布斯采样
玻尔兹曼机的吉布斯采样过程为：随机选择一个变量Xi，然后根据其全条件概率P(Xi|X-i)来设置其状态值，即以P(Xi=1|X-i)的概率将变量Xi设为1，否则为0。在固定的温度T下，运行足够时间后，玻尔兹曼机会达到热平衡状态。此时，任何全局状态的概率都服从玻尔兹曼分布P(x)，只和系统的能量有关，和初始状态无关。

如果你听不懂上述定义，(针对没有概率论和热力学基础的读者)。那你只需要了解，因为玻尔兹曼机的联合概率函数中配分函数Z很难处理，就采用另一种方法来将节点的状态值概率趋近于玻尔兹曼分布。

玻尔兹曼机可以解决两类问题，一类是搜索问题：当给定变量之间的连接权重时，需要找到一组二值向量，使得整个网络的能量最低。另一类是学习问题，当给定变量的多组观测值时，学习网络的最优权重。

受限玻尔兹曼机(RBM)

全连接的玻尔兹曼机在理论上固然有趣，但是由于其复杂性，目前为止并没有广泛运用。实际应用中，用得更多的是基于玻尔兹曼机改造的一个版本——受限玻尔兹曼机(RBM)，其网络架构如下：
在这里插入图片描述
玻尔兹曼机没有层的概念，它所有的节点都是全连接的。
但受限玻尔兹曼机有层的概念。它有两层，一层称为显层，用于观测和输入，一层为隐藏层，用于提取特征。
受限玻尔兹曼机相比玻尔兹曼机，层间的节点还是采用了对称的全连接的方式连接，但是层内的节点相互独立，相互不受影响。
因为层内节点相互独立，那么由Bayes条件独立定理，受限玻尔兹曼机可以并行地对所有的显层变量或隐藏层变量同时进行采样，从而更快达到热平衡。

能量函数

在这里插入图片描述
由于受限玻尔兹曼机变成了层的结构，所以受限玻尔兹曼机的能量函数也变成了由三部分组成，
一个是显层节点偏置乘以显层随机可观测变量部分，
一个是连接权重与显层随机可观测变量和隐层随机可观测变量相乘部分，
一个是隐层节点偏置乘以隐层随机可观测变量偏置部分。

CD学习算法

由于受限玻尔兹曼机的特殊结构，G·Hinton提出了一种比吉布斯采样更加有效的学习算法，即对比散度学习算法，又称为CD学习算法。
通过对CD学习算法的学习，我发现这个CD算法就是在吉布斯采样的基础上作出的一点改进，即在处理玻尔兹曼机时，运行无穷次的吉布斯采样改进为运行K次即可。
以前处理玻尔兹曼机时，吉布斯采样是一直对这个玻尔兹曼机处理，直到这个波尔兹曼机收敛。G·Hinton提出，在受限玻尔兹曼机中，不需要等到受限玻尔兹曼机完全收敛，只需要K步吉布斯采样，这时模型就非常好了。所以CD算法又称K步吉布斯采样法。

代码实现受限玻尔兹曼机

受限玻尔兹曼机有两个偏置项，隐藏层的偏置项有助于RBM在前向传递中获得非零激活值，而可见层的偏置项有助于受限玻尔兹曼机学习后向传递中的重建。
在正向传递中，每个输入数据乘以一个独立的权重，然后相加后再加上一个偏置项，最后将结果传递到激活函数来产生输出。
用对比散度计算正反向梯度，然后更新偏置和权重
因为最开始受限玻尔兹曼机权重是随机初始化的，所以重建结果和原始输入差距通常会比较大，这个差距可看作是重建误差，训练受限玻尔兹曼机是通过在可见层和隐藏层之间迭代学习不断正向反向传播，直至达到某个误差的最小值
我们采用了Python环境下Numpy工具库来撰写代码，具体注释已在代码中标注，不做过多讲解。

import numpy
class RBM:
    def __init__(self, n_visible, n_hidden):
        self.n_visible = n_visible
        self.n_hidden = n_hidden
        self.bias_a = np.zeros(self.n_visible)  # 可视层偏移量
        self.bias_b = np.zeros(self.n_hidden)  # 隐藏层偏移量
        self.weights = np.random.normal(0, 0.01, size=(self.n_visible, self.n_hidden))
        self.n_sample = None

    def encode(self, v):
        # 编码，即基于v计算h的条件概率：p(h=1|v)
        return sigmoid(self.bias_b + v @ self.weights)

    def decode(self, h):
        # 解码(重构)：即基于h计算v的条件概率：p(v=1|h)
        return sigmoid(self.bias_a + h @ self.weights.T)

    def gibbs_sample(self, v0, max_cd):
        # gibbs采样, 返回max_cd采样后的v以及h值
        v = v0
        for _ in range(max_cd):
            # 首先根据输入样本对每个隐藏层神经元采样。二项分布采样，决定神经元是否激活
            ph = self.encode(v)
            h = np.random.binomial(1, ph, (self.n_sample, self.n_hidden))
            # 根据采样后隐藏层神经元取值对每个可视层神经元采样
            pv = self.decode(h)
            v = np.random.binomial(1, pv, (self.n_sample, self.n_visible))
        return v
        
    def update(self, v0, v_cd, eta):
        # 根据Gibbs采样得到的可视层取值(解码或重构)，更新参数
        ph = self.encode(v0)
        ph_cd = self.encode(v_cd)
        self.weights += eta * (v0.T @ ph - v_cd.T @ ph)  # 更新连接权重参数
        self.bias_b += eta * np.mean(ph - ph_cd, axis=0)  # 更新隐藏层偏移量b
        self.bias_a += eta * np.mean(v0 - v_cd, axis=0)  # 更新可视层偏移量a
        return

    def train(self, data, max_step, max_cd=2, eta=0.1):
        # 训练主函数,采用对比散度算法(CD算法)更新参数
       
        assert data.shape[1] == self.n_visible, "输入数据维度与可视层神经元数目不相等"
        self.n_sample = data.shape[0]
        for i in range(max_step):
            v_cd = self.gibbs_sample(data, max_cd)
            self.update(data, v_cd, eta)
            error = np.sum((data - v_cd) ** 2) / self.n_sample / self.n_visible * 100
            if i == (max_step-1):  # 将重构后的样本与原始样本对比计算误差
                print("可视层(隐藏层)状态误差比例:{0}%".format(round(error, 2)))
    def predict(self, v):
        # 输入训练数据，预测隐藏层输出
        ph = self.encode(v)[0]
        states = ph >= np.random.rand(len(ph))
        return states.astype(int)

深度玻尔兹曼机(DBM)

在这里插入图片描述

深度玻尔兹曼机实际上是由多个受限玻尔兹曼机堆栈构成，我们构建一个简单的三层深度玻尔兹曼机，就首先需要进行预训练，通过对比散度算法训练出两个受限玻尔兹曼机，对于每个受限玻尔兹曼机，都是根据可见层的数据，来学习到隐层的数据，对于第二个受限玻尔兹曼机，可以把求得的第一个隐层看作是可见层来求第二个隐层的数据，然后再combine两个受限玻尔兹曼机进行微调，通过CD算法来更新相应的数据，最后实现一个简单的三层深度玻尔兹曼机。
对于受限玻尔兹曼机每新增的隐藏层，权重都会通过迭代学习反复调整，直至该层能够逼近前一层的输入，这是贪婪的、逐层的无监督的预训练。
深度玻尔兹曼机网络结构如下：
在这里插入图片描述

代码实现深度玻尔兹曼机

# Implementation of 3 layer ( 2 hidden layer ) Deep Boltzmann Machine
import numpy as np     
from scipy.special import expit   #explit函数也称为logistic sigmoid函数，expit（x）= 1 /（1 + exp（-x））
from matplotlib import pylab as plt    

def binary_cross_entropy(data, reconst):    #交叉熵损失函数(交叉熵就是用来判定实际的输出与期望的输出的接近程度)
    return - np.mean( np.sum( data * np.log(reconst) + (1-data) * np.log(1 - reconst), axis=1) )

def reconstruct_data(data, b, c1, c2, w_vh1, w_h1h2, num_sample=100): #重建数据
    m_h1 = expit( np.dot(data, w_vh1) + c1 )       #自下而上求得第一个隐层的数据，dot()函数用来做矩阵乘法
    for i in range(num_sample):
        m_h2 = expit( np.dot(m_h1, w_h1h2) + c2 )  #自下而上求得第二个隐层的数据
        m_h1 = expit( np.dot(w_h1h2, m_h2.T).T + c1 )  #第一个隐层的重建数据
    return expit( np.dot(w_vh1, m_h1.T).T + b )     #返回可视层的重建数据

def popup(data, c1, w_vh1):     #编码函数
    return expit(np.dot(data, w_vh1) + c1)

def rbm_contrastive_divergence(data, b, c, w, num_sample=100):  #RBM的对比散度算法(CD算法)
    # Mean field可视层参数和隐藏层参数定义
    m_vis = data
    m_hid = expit( np.dot(data, w) + c )
     
    s_vis = m_vis
    for i in range(num_sample):
        #Gibbs采样，决定神经元是否被激活
        sm_hid = expit( np.dot(s_vis, w) + c )
        s_hid = np.random.binomial(1, sm_hid)       #吉布斯采样后的隐藏层数据
        sm_vis = expit( np.dot(w, s_hid.T).T + b )
        s_vis = np.random.binomial(1, sm_vis)	#吉布斯采样后的可见层数据
    return np.mean(m_vis - s_vis, axis=0), np.mean(m_hid - s_hid, axis=0), 
                    (np.dot(m_vis.T, m_hid) - np.dot(s_vis.T, s_hid)) / len(data)    #返回原始数据和最后一次采样数据梯度的平均值

def dbm_contrastive_divergence(data, b, c1, c2, w_vh1, w_h1h2, num_sample=100):		#相比于RBM多一个隐藏层的DBM的对比散度算法(CD算法)
    # Mean field相应参数的定义
    m_vis = data
    m_h1 = np.random.uniform(size=(len(data), len(c1)))    
    m_h2 = np.random.uniform(size=(len(data), len(c2)))
    for i in range(num_sample):
        m_h1 = expit( np.dot(m_vis, w_vh1) + np.dot(w_h1h2, m_h2.T).T + c1 )
        m_h2 = expit( np.dot(m_h1, w_h1h2) + c2 )
    # Gibbs sample吉布斯采样
    s_vis = np.random.binomial(1, m_vis)
    s_h1 = np.random.binomial(1, 0.5, size=(len(data), len(c1)))
    s_h2 = np.random.binomial(1, 0.5, size=(len(data), len(c2)))
    for i in range(num_sample):
        sm_vis = expit( np.dot(w_vh1, s_h1.T).T + b )
        s_vis = np.random.binomial(1, sm_vis)
        sm_h1 = expit( np.dot(s_vis, w_vh1) + np.dot(w_h1h2, s_h2.T).T + c1 )
        s_h1 = np.random.binomial(1, sm_h1)
        sm_h2 = expit( np.dot(s_h1, w_h1h2) + c2 )
        s_h2 = np.random.binomial(1, sm_h2)
    return np.mean(m_vis - s_vis, axis=0), np.mean(m_h1 - s_h1, axis=0), np.mean(m_h2 - s_h2, axis=0), 
                ( np.dot(m_vis.T, m_h1) - np.dot(s_vis.T, s_h1) ) / len(data), ( np.dot(m_h1.T, m_h2) - np.dot(s_h1.T, s_h2) ) / len(data)

# Assign structural parameters 分配结构参数
num_visible = 784   #可视层神经元节点个数
num_hidden1 = 500   #第一个隐藏层神经元节点个数
num_hidden2 = 1000  #第二个隐藏层神经元节点个数

# Assign learning parameters分配学习参数
pretrain_epochs = 100           #预训练的迭代次数
pretrain_learning_rate = 0.1    #预训练学习率
train_epochs = 100              #训练迭代次数
train_learning_rate = 0.1       #训练学习率

# Initialize weights and biases  初始化权重，偏置
b = np.zeros((num_visible, ))      
c1 = np.zeros((num_hidden1, ))
c2 = np.zeros((num_hidden2, ))
 #随机初始化权重
w_vh1 = np.random.normal(scale=0.01, size=(num_visible, num_hidden1))     
w_h1h2 = np.random.normal(scale=0.01, size=(num_hidden1, num_hidden2))

# Load data, data needs to be in range [0, 1]   
data = np.load("./imgs.npy").reshape(10, 28*28)[[3, 7, 9]]

# Pretraining预训练
for i in range(pretrain_epochs):
    # Calculate gradient 计算梯度
    update_b, update_c1, update_w_vh1 = rbm_contrastive_divergence(data, b, c1, w_vh1)
    # Upate parameters 更新可视层和第一个隐层的偏置以及连接权
    b += pretrain_learning_rate * update_b
    c1 += pretrain_learning_rate * update_c1
    w_vh1 += pretrain_learning_rate * update_w_vh1

pseudo_data = popup(data, c1, w_vh1) #保存由可见层编码学习到的第一个隐层的数据
for i in range(pretrain_epochs):
    # Calculate gradient计算梯度
    update_c1, update_c2, update_w_h1h2 = rbm_contrastive_divergence(pseudo_data, c1, c2, w_h1h2)
    # Upate parameters更新两个隐层的偏置以及连接权
    c1 += pretrain_learning_rate * update_c1
    c2 += pretrain_learning_rate * update_c2
    w_h1h2 += pretrain_learning_rate * update_w_h1h2

# Show current cost 打印原始数据与重建结果之间的误差
cost = binary_cross_entropy(data, reconstruct_data(data, b, c1, c2, w_vh1, w_h1h2))
print( "Reconstruction cost is %.2f"%cost )

# Fine tuning 微调
for i in range(train_epochs):
    # Calculate gradient
    update_b, update_c1, update_c2, update_w_vh1, update_w_h1h2 
                                    = dbm_contrastive_divergence(data, b, c1, c2, w_vh1, w_h1h2)
    # Update parameters更新所有的偏置以及连接权
    b += train_learning_rate * update_b
    c1 += train_learning_rate * update_c1
    c2 += train_learning_rate * update_c2
    w_vh1 += train_learning_rate * update_w_vh1
    w_h1h2 += train_learning_rate * update_w_h1h2

# Show fine tuning result 打印微调后的原始数据与重建结果之间的误差
cost = binary_cross_entropy(data, reconstruct_data(data, b, c1, c2, w_vh1, w_h1h2))
print( "Reconstruction cost is %.2f"%cost )

# Show result images 显示.npy文件最后形成的灰度图
plt.matshow(reconstruct_data(data, b, c1, c2, w_vh1, w_h1h2)[0].reshape(28, 28))
plt.gray()
plt.show()

深度置信网络(DBN)

玻尔兹曼机和深度置信网络是生成模型，借助隐变量来描述复杂的数据分布。
深度置信网络是神经网络的一种。既可以用于非监督学习，也可以用于监督学习。
深度置信网络组成元件是受限玻尔兹曼机。通过下图的网络结构，我们可以看出深度置信网络和受限玻尔兹曼机的关系：
在这里插入图片描述
训练深度置信网络由两部分组成。
一是单独训练每一个受限玻尔兹曼机，使它们收敛。
二是将每一个受限玻尔兹曼机深入展开，构成一个前向传播的深层网络。然后使用BP算法微调参数，从而使整个网络收敛。

代码实现深度置信网络

我们采用了Python环境下Numpy工具库来和tensorflow2来撰写代码，具体注释已在代码中标注，不做过多讲解。

import numpy as np
import tensorflow as tf

def sigmoid(z):
    return 1 / (1 + np.exp(-z))

def DBN(epoch, test_data, rbm_list,  lr=0.001):
    # 构建深度置信网络

    RBM1_w = tf.Variable(tf.convert_to_tensor(rbm_list[0].weights))     # 将数据转换为Varible类型，方便梯度跟踪
    RBM1_vb = tf.Variable(tf.convert_to_tensor(rbm_list[0].bias_a))     # 将数据转换为Varible类型，方便梯度跟踪
    RBM1_hb = tf.Variable(tf.convert_to_tensor(rbm_list[0].bias_b))     # 将数据转换为Varible类型，方便梯度跟踪

    RBM1to2_w = tf.Variable(tf.random.normal([6, 6],dtype=tf.float64))  # 将数据转换为Varible类型，方便梯度跟踪

    RMB2_w = tf.Variable(tf.convert_to_tensor(rbm_list[1].weights))     # 将数据转换为Varible类型，方便梯度跟踪
    RBM2_vb = tf.Variable(tf.convert_to_tensor(rbm_list[1].bias_a))     # 将数据转换为Varible类型，方便梯度跟踪
    RBM2_hb = tf.Variable(tf.convert_to_tensor(rbm_list[1].bias_b))     # 将数据转换为Varible类型，方便梯度跟踪

    BP_w = tf.Variable(tf.random.normal([6,6],dtype=tf.float64))        # 将数据转换为Varible类型，方便梯度跟踪
    BP_b = tf.Variable(tf.random.normal([6],dtype=tf.float64))          # 将数据转换为Varible类型，方便梯度跟踪

    test_data = tf.Variable(tf.convert_to_tensor(test_data,dtype=tf.float64))

    # 微调迭代次数
    for step in range(epoch):
        # 计算当前批次样本的网络前向传播
        with tf.GradientTape() as tape:
            # 前向传播
            # 第一层计算 第一个RBM v-h
            out1 = tf.nn.sigmoid(RBM1_hb + test_data @ RBM1_w)
            # 第二层计算 h-v
            out2 = tf.nn.sigmoid(RBM2_vb + out1 @ RBM1to2_w)
            # 第三层输出也就是第二个RBM v-h
            out3 = tf.nn.sigmoid(RBM2_hb + out2 @ RMB2_w)
            out = tf.nn.relu(BP_b + out3 @ BP_w)
            # 反向传播
            # 计算损失函数
            loss = tf.reduce_mean(tf.square(test_data - out))   # 均方差损失函数
            # 手动梯度更新参数
            grads = tape.gradient(loss, [RBM1_w, RBM1_vb, RBM1_hb,
                                         RBM1to2_w,
                                         RMB2_w, RBM2_vb, RBM2_hb,
                                         BP_b,BP_w])
            # 参数更新
            RBM1_w.assign_sub(lr * grads[0])
            # RBM1_vb = RBM1_vb - lr * grads[1]
            RBM1_hb.assign_sub(lr * grads[2])
            RBM1to2_w.assign_sub(lr * grads[3])
            RMB2_w.assign_sub(lr * grads[4])
            RBM2_vb.assign_sub(lr * grads[5])
            RBM2_hb.assign_sub(lr * grads[6])
            BP_b.assign_sub(lr * grads[7])
            BP_w.assign_sub(lr * grads[8])

            print(step,": loss : ", loss)

class RBM:
    def __init__(self, n_visible, n_hidden):
        self.n_visible = n_visible          # 可视层节点个数
        self.n_hidden = n_hidden            # 隐藏层节点个数
        self.bias_a = np.zeros(self.n_visible)  # 可视层偏移量
        self.bias_b = np.zeros(self.n_hidden)  # 隐藏层偏移量
        self.weights = np.random.normal(0, 0.01, size=(self.n_visible, self.n_hidden))  # 连接权重w
        self.n_sample = None

    def encode(self, v):
        # 编码，即基于v计算h的条件概率：p(h=1|v)
        return sigmoid(self.bias_b + v @ self.weights)

    def decode(self, h):
        # 解码(重构)：即基于h计算v的条件概率：p(v=1|h)
        return sigmoid(self.bias_a + h @ self.weights.T)

    # gibbs采样, 返回max_cd采样后的v以及h值
    def gibbs_sample(self, v0, max_cd):
        v = v0
        for _ in range(max_cd):
            # 首先根据输入样本对每个隐藏层神经元采样。二项分布采样，决定神经元是否激活
            ph = self.encode(v)
            h = np.random.binomial(1, ph, (self.n_sample, self.n_hidden))
            # 根据采样后隐藏层神经元取值对每个可视层神经元采样
            pv = self.decode(h)
            v = np.random.binomial(1, pv, (self.n_sample, self.n_visible))
        return v

    # 根据Gibbs采样得到的可视层取值(解码或重构)，更新参数
    def update(self, v0, v_cd, eta):
        ph = self.encode(v0)
        ph_cd = self.encode(v_cd)
        self.weights += eta * (v0.T @ ph - v_cd.T @ ph)  # 更新连接权重参数
        self.bias_b += eta * np.mean(ph - ph_cd, axis=0)  # 更新隐藏层偏移量b
        self.bias_a += eta * np.mean(v0 - v_cd, axis=0)  # 更新可视层偏移量a
        return

    # 训练函数 采用对比散度算法更新参数
    def fit(self, data, max_step, max_cd=2, eta=0.1):
        # data 训练数据集
        # max_cd 采样步数
        # max_step: 最大迭代次数 iter
        # eta: 学习率
        assert data.shape[1] == self.n_visible, "输入数据维度与可视层神经元数目不相等"
        self.n_sample = data.shape[0]

        for i in range(max_step):
            v_cd = self.gibbs_sample(data, max_cd)
            self.update(data, v_cd, eta)
            error = np.sum((data - v_cd) ** 2) / self.n_sample / self.n_visible * 100
            if i == (max_step-1):  # 将重构后的样本与原始样本对比计算误差
                print("可视层(隐藏层)状态误差比例:{0}%".format(round(error, 2)))

    # 预测
    def predict(self, v):
        # 输入训练数据，预测隐藏层输出
        ph = self.encode(v)[0]
        states = ph >= np.random.rand(len(ph))
        return states.astype(int)


if __name__ == '__main__':
    # 迭代次数
    iter = 100
    # 学习率
    lr = 0.001
    # N代表组成DBN的RBM的层数
    N = 2

    # 创建多个RBM层
    rbm_model_list = []
    # 用于两个RBM的连接权重，h-v
    for i in range(N):
        rbm_model_list.append(RBM(n_visible=6, n_hidden=6))

    # 训练集数据
    V = np.array([[1, 1, 1, 0, 0, 0], [1, 0, 1, 0, 0, 0], [1, 1, 1, 0, 0, 0],
                  [0, 0, 1, 1, 1, 0], [0, 0, 1, 1, 0, 0], [0, 0, 1, 1, 1, 0],
                  [1, 0, 0, 1, 1, 0], [0, 0, 0, 1, 0, 0], [0, 0, 1, 0, 1, 0]])

    # 单独训练每一个RBM层的参数
    for epoch in range(iter):   # 迭代次数
        for i in range(N):      # 每一次迭代 单独训练每一个RBM层
            rbm_model_list[i].fit(V, max_step=iter, max_cd=1, eta=lr)

    # 将两层RBM串联起来进行预测
    user = np.array([[1,0,0,1,0,0]])
    temp = rbm_model_list[0].predict(user)
    out = rbm_model_list[1].predict([temp])
    print(out)

    DBN(iter,V,rbm_model_list)