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)