Question of about Boltzmann machine implementation

Question

I would like to implement Boltzmann machine with two hidden and two visible units. The four possible hidden units configuration are $(0,0), (0,1), (1,0), (1,1)$, and their probability distribution is $20\%, 35\%, 35\%, 10\%$, respectively. I would like to train Boltzmann machine so that it learns this distribution. The architecture is there are connections between hidden-hidden, visible-hidden units but no connection exists between two visible units. Here is a code that seems to be designed for this problem, the thing that I am confused about is how do I feed the probability distribution data to the training, and how do I test if the algorithm is trained well for this distribution.

#!/usr/bin/env python
from future import division
from enum import Enum
import numpy as np
Clamp = Enum('Clamp', 'VISIBLE_UNITS NONE INPUT_UNITS')
class Step:
    def init(self, temperature, epochs):
        self.temperature = temperature
        self.epochs = epochs
numInputUnits = 2
numHiddenUnits = 2
numVisibleUnits = numInputUnits 
numUnits = numVisibleUnits+numHiddenUnits
annealingSchedule = [Step(20.,2),
                     Step(15.,2),
                     Step(12.,2),
                     Step(10.,4)]
coocurranceCycle = Step(10.,10)
weights = np.zeros((numUnits,numUnits))
states = np.zeros(numUnits)
energy = np.zeros(numUnits)
connections = np.zeros((numUnits,numUnits), dtype=np.int)
for i in range(numInputUnits):
    for j in range(1,numHiddenUnits+1):
        connections[i,-j] = 1
for i in range(numHiddenUnits,0,-1):
    for j in range(i-1,0,-1):
        connections[-i,-j] = 1
valid = np.nonzero(connections)
numConnections = np.size(valid[0])
connections[valid] = np.arange(1,numConnections+1)
connections = connections + connections.T - 1
def propagate(temperature, clamp):
    global energy, states, weights
if clamp == Clamp.VISIBLE_UNITS:
    numUnitsToSelect = numHiddenUnits
elif clamp == Clamp.NONE:
    numUnitsToSelect = numUnits
else:
    numUnitsToSelect = numHiddenUnits

for i in range(numUnitsToSelect):
    # Calculating the energy of a randomly selected unit    
    unit=numUnits-np.random.randint(1,numUnitsToSelect+1)
    energy[unit] = np.dot(weights[unit,:], states)

    p = 1. / (1.+ np.exp(-energy[unit] / temperature))
    states[unit] = 1. if  np.random.uniform() <= p else 0 


def anneal(annealingSchedule, clamp):
    for step in annealingSchedule:
        for epoch in range(step.epochs):
            propagate(step.temperature, clamp)
def sumCoocurrance(clamp):

    sums = np.zeros(numConnections)
    for epoch in range(coocurranceCycle.epochs):
        propagate(coocurranceCycle.temperature, clamp)
        for i in range(numUnits):
            if(states[i] == 1):
                for j in range(i+1,numUnits):
                    if(connections[i,j]>-1 and states[j] ==1):
                        sums[connections[i,j]] += 1

    return sums
def updateWeights(pplus, pminus):
    global weights
    for i in range(numUnits):
        for j in range(i+1,numUnits):

            if connections[i,j] > -1:
                index = connections[i,j]
                weights[i,j] += 2*np.sign(pplus[index] - pminus[index])
                weights[j,i] = weights[i,j]
def recall(pattern):
    global states
# Setting pattern to recall
states[0:numInputUnits] = pattern

# Assigning random values to the hidden and output states
states[-(numHiddenUnits):] = np.random.choice([0,1],numHiddenUnits)

anneal(annealingSchedule, Clamp.INPUT_UNITS)

return states[numInputUnits:numInputUnits]


def addNoise(pattern):
    probabilities = 0.8*pattern+0.05
    uniform = np.random.random(numVisibleUnits)

    return (uniform < probabilities).astype(int)
def learn(patterns):
    global states, weights
numPatterns = patterns.shape[0]    
trials=numPatterns*coocurranceCycle.epochs
weights = np.zeros((numUnits,numUnits))

for i in range(1800):
    # Positive phase
    pplus = np.zeros(numConnections)
    for pattern in patterns:

        # Setting visible units values (inputs and outputs)
        states[0:numVisibleUnits] = addNoise(pattern)

        # Assigning random values to the hidden units
        states[-numHiddenUnits:] = np.random.choice([0,1],numHiddenUnits)

        anneal(annealingSchedule, Clamp.VISIBLE_UNITS)
        pplus += sumCoocurrance(Clamp.VISIBLE_UNITS)
    pplus /= trials

    # Negative phase
    states = np.random.choice([0,1],numUnits)          
    anneal(annealingSchedule, Clamp.NONE)
    pminus = sumCoocurrance(Clamp.NONE) / coocurranceCycle.epochs

    updateWeights(pplus,pminus)