diff --git a/demos/CMakeLists.txt b/demos/CMakeLists.txt
index fcc7d0b8..98f06e3a 100644
--- a/demos/CMakeLists.txt
+++ b/demos/CMakeLists.txt
@@ -72,3 +72,6 @@ create_demo(Demo_Stereo "Demo Stereo" "DGM")
 
 # ================================================ DEMO ICR =================================================
 create_demo(Demo_ICR "Demo ICR" "DGM;DNN;VIS")
+
+# ================================================ DEMO RBM =================================================
+create_demo(Demo_RBM "Demo RBM" "DGM;DNN;VIS")
diff --git a/demos/Demo RBM.cpp b/demos/Demo RBM.cpp
new file mode 100644
index 00000000..bba45315
--- /dev/null
+++ b/demos/Demo RBM.cpp	
@@ -0,0 +1,124 @@
+#include "DNN.h"
+#include "DGM.h"
+#include "VIS.h"
+#include "DGM/timer.h"
+#include <fstream>
+
+namespace dgm = DirectGraphicalModels;
+
+/**
+ * Reads the digits numerical value in a decimal notation
+ *
+ * @param file to read, and the number of digits to read
+ * @return an array of digit labels
+ */
+std::vector<byte> readGroundTruth(const std::string& fileName)
+{
+	std::vector<byte> res;
+	std::ifstream inFile;
+	inFile.open(fileName.c_str());
+
+	if (inFile.is_open()) {
+		int val;
+		while (!inFile.eof()) {
+			inFile >> val;
+			res.push_back(static_cast<byte>(val));
+		}
+		inFile.close();
+	}
+	return res;
+}
+
+
+float sigmoidFunction(float x)
+{
+	return 1.0f / (1.0f + expf(-x));
+}
+
+float sigmoidFunction_derivative(float x)
+{
+	float s = sigmoidFunction(x);
+	return s * (1 - s);
+}
+
+int main()
+{
+	const float learningRate             = 0.05f;
+	const size_t numEpochs				 = 20;
+	const size_t numTestSamples			 = 2000;
+	const size_t numTrainSamples		 = 4000;
+	
+	//const byte	nStates					 = 10;
+	const word nFeatures                 = 28 * 28;
+	const size_t numNeuronsHiddenLayer   = 10;
+
+#ifdef WIN32
+	const std::string dataPath = "../../data/digits/";
+#else
+	const std::string dataPath = "../../../data/digits/";
+#endif
+
+	auto pLayerVisible = std::make_shared<dgm::dnn::CNeuronLayer>(nFeatures, 1, [](float x) { return x; }, [](float x) { return 1.0f; });
+	auto pLayerHidden  = std::make_shared<dgm::dnn::CNeuronLayer>(numNeuronsHiddenLayer, nFeatures, &sigmoidFunction, &sigmoidFunction_derivative);
+
+	pLayerVisible->generateRandomWeights();
+	pLayerHidden->generateRandomWeights();
+	
+	dgm::dnn::CRBM rbm({ pLayerVisible, pLayerHidden });
+
+	//rbm.debug();
+	Mat fv;
+
+	// ==================== TRAINING DIGITS ====================
+	dgm::Timer::start("Training...\n");
+	auto	trainGT = readGroundTruth(dataPath + "train_gt.txt");
+	for (size_t e = 0; e < numEpochs; e++)
+		for (int s = 0; s < numTrainSamples; s++) {
+			std::stringstream ss;
+			ss << dataPath << "train/digit_" << std::setfill('0') << std::setw(4) << s << ".png";
+			std::string fileName = samples::findFile(ss.str());
+			Mat img = imread(fileName, 0);
+			img = img.reshape(1, img.cols * img.rows);
+			img.convertTo(fv, CV_32FC1, 1.0 / 255);
+			fv = Scalar(1.0f) - fv;
+
+			rbm.contrastiveDivergence(fv, 0.5f);
+		} // samples
+	dgm::Timer::stop();
+
+	// ==================== TESTING DIGITS ====================
+	//dgm::CCMat confMat(nStates);
+	dgm::Timer::start("Testing...");
+	auto 	testGT = readGroundTruth(dataPath + "test_gt.txt");
+	for (size_t s = 0; s < numTestSamples; s++) {
+		std::stringstream ss;
+		ss << dataPath << "test/digit_" << std::setfill('0') << std::setw(4) << s << ".png";
+		std::string fileName = samples::findFile(ss.str());
+		Mat img = imread(fileName, 0);
+		img = img.reshape(1, img.cols * img.rows);
+		img.convertTo(fv, CV_32FC1, 1.0 / 255);
+		fv = Scalar(1.0f) - fv;
+
+		Mat outputValues = rbm.reconstruct(fv);
+
+		//Point maxclass;
+		//minMaxLoc(outputValues, NULL, NULL, NULL, &maxclass);
+		//int number = maxclass.y;
+
+		//confMat.estimate(number, testGT[s]);
+		//printf("prediction [%d] for digit %d with %.3f%s at position %zu \n", number, testDataDigit[z], maxAccuracy, "%", z);
+	} // samples
+	dgm::Timer::stop();
+	//printf("Accuracy = %.2f%%\n", confMat.getAccuracy());
+
+	// Confusion matrix
+	//dgm::vis::CMarker marker;
+	//Mat cMat = confMat.getConfusionMatrix();
+	//Mat cMatImg = marker.drawConfusionMatrix(cMat, dgm::vis::MARK_BW);
+	//imshow("Confusion Matrix", cMatImg);
+	rbm.debug();
+
+	waitKey();
+
+	return 0;
+}
diff --git a/include/DNN.h b/include/DNN.h
index 7eb05db3..7051d58a 100644
--- a/include/DNN.h
+++ b/include/DNN.h
@@ -2,6 +2,7 @@
 #include "DNN/Neuron.h"
 #include "DNN/NeuronLayer.h"
 #include "DNN/Perceptron.h"
+#include "DNN/RBM.h"
 // #include "DNN/Functions.hpp"
 
 /**
diff --git a/modules/DNN/CMakeLists.txt b/modules/DNN/CMakeLists.txt
index 6e8924de..bff0659c 100644
--- a/modules/DNN/CMakeLists.txt
+++ b/modules/DNN/CMakeLists.txt
@@ -9,6 +9,7 @@ source_group("" FILES ${DNN_SOURCES} ${DNN_HEADERS})
 source_group("Source Files\\Neuron" FILES	"Neuron.h" "Neuron.cpp")
 source_group("Source Files\\Neuron Layer" FILES	"NeuronLayer.h" "NeuronLayer.cpp")
 source_group("Source Files\\Perceptron" FILES	"Perceptron.h" "Perceptron.cpp")
+source_group("Source Files\\RBM" FILES	"RBM.h" "RBM.cpp")
 
 # Properties -> C/C++ -> General -> Additional Include Directories
 include_directories(${PROJECT_SOURCE_DIR}/include
diff --git a/modules/DNN/NeuronLayer.cpp b/modules/DNN/NeuronLayer.cpp
index 4009c0b6..60c65160 100644
--- a/modules/DNN/NeuronLayer.cpp
+++ b/modules/DNN/NeuronLayer.cpp
@@ -1,39 +1,45 @@
-#include "NeuronLayer.h"
-#include "DGM/random.h"
-#include "DGM/parallel.h"
-#include "macroses.h"
-
-namespace DirectGraphicalModels { namespace dnn
-{
-	void CNeuronLayer::generateRandomWeights(void)
-	{
-		m_weights = random::U(m_weights.size(), m_weights.type(), -0.5f, 0.5f);
-		m_biases = random::U(m_biases.size(), m_biases.type(), -0.5f, 0.5f);
-	}
-
-	void CNeuronLayer::dotProd(const Mat& values)
-	{
-		// this->m_netValues = this->m_weights * values + m_biases;
-		gemm(m_weights.t(), values, 1, m_biases, 1, m_netValues);
-	}
-
-	void CNeuronLayer::setNetValues(const Mat& values)
-	{
-		// Assertions
-		DGM_ASSERT(values.type() == m_netValues.type());
-		DGM_ASSERT(values.size() == m_netValues.size());
-		values.copyTo(m_netValues);
-	}
-
-	Mat CNeuronLayer::getValues(void) const
-	{ 
-		Mat res(m_netValues.clone());
-		for (int y = 0; y < res.rows; y++) {
-			float* pRes = res.ptr<float>(y);
-			for (int x = 0; x < res.cols; x++)
-				pRes[x] = m_activationFunction(pRes[x]);
-		}
-		return res; 
-	}
-
-}}
+#include "NeuronLayer.h"
+#include "DGM/random.h"
+#include "DGM/parallel.h"
+#include "macroses.h"
+
+namespace DirectGraphicalModels {
+	namespace dnn
+	{
+		void CNeuronLayer::generateRandomWeights(void)
+		{
+			m_weights = random::U(m_weights.size(), m_weights.type(), -0.5f, 0.5f);
+			m_biases = random::U(m_biases.size(), m_biases.type(), -0.5f, 0.5f);
+		}
+
+		void CNeuronLayer::dotProd(const Mat& values)
+		{
+			// this->m_netValues = this->m_weights * values + m_biases;
+			gemm(m_weights.t(), values, 1, m_biases, 1, m_netValues);
+		}
+		void CNeuronLayer::dotProdVis(const Mat& values, const Mat& weights)
+		{
+			// this->m_netValues = weights * values + m_biases;
+			gemm(weights, values, 1, m_biases, 1, m_netValues);
+		}
+		void CNeuronLayer::setNetValues(const Mat& values)
+		{
+			// Assertions
+			DGM_ASSERT(values.type() == m_netValues.type());
+			DGM_ASSERT(values.size() == m_netValues.size());
+			values.copyTo(m_netValues);
+		}
+
+		Mat CNeuronLayer::getValues(void) const
+		{
+			Mat res(m_netValues.clone());
+			for (int y = 0; y < res.rows; y++) {
+				float* pRes = res.ptr<float>(y);
+				for (int x = 0; x < res.cols; x++)
+					pRes[x] = m_activationFunction(pRes[x]);
+			}
+			return res;
+		}
+
+	}
+}
\ No newline at end of file
diff --git a/modules/DNN/NeuronLayer.h b/modules/DNN/NeuronLayer.h
index 72de381e..7b7b2fec 100644
--- a/modules/DNN/NeuronLayer.h
+++ b/modules/DNN/NeuronLayer.h
@@ -38,6 +38,13 @@ namespace DirectGraphicalModels {
 			*/
 			DllExport void  dotProd(const Mat& values);
 			/**
+			* In the dotProd method, we cant multiply hidden neuron values by hidden neuron weights, so dotProdVis is created
+			* which we can input by which weights neuron are multiplied by. 
+			* in dotProd ->     this->m_netValues = this->m_weights * values + this->m_biases;
+			* in dotProdVis ->  this->m_netValues = m_weights * values + this->m_biases;
+			*/
+			DllExport void  dotProdVis(const Mat& values, const Mat& weights); 
+			/**
 			* @brief Returns the values of the neurons of the layer
 			* @note This method returns the result of per-element application of the activation function to the neurons' net values, i.e. activationFunction(netValues)
 			* @returns The values of the neurons of the layer (size: 1 x numNeurons; type: CV_32FC1)
@@ -64,4 +71,3 @@ namespace DirectGraphicalModels {
 		using ptr_nl_t = std::shared_ptr<CNeuronLayer>;
 	}
 }
-
diff --git a/modules/DNN/RBM.cpp b/modules/DNN/RBM.cpp
new file mode 100644
index 00000000..7f6a05f3
--- /dev/null
+++ b/modules/DNN/RBM.cpp
@@ -0,0 +1,121 @@
+#include "RBM.h"
+#include "DGM/random.h"
+#include "macroses.h"
+
+namespace DirectGraphicalModels {
+	namespace dnn {
+		CRBM::CRBM(const std::vector<ptr_nl_t>& vpLayers){
+			for (auto& nl : vpLayers)
+				m_vpNeuronLayers.push_back(nl);
+		}
+
+		Mat CRBM::getBinomial(const Mat& mean) {
+			Mat res(mean.clone());
+			for (int y = 0; y < res.rows; y++) {
+				float* pRes = res.ptr<float>(y);
+				for (int x = 0; x < res.cols; x++) {
+					if (pRes[x] < 0 || pRes[x]>1) {
+						pRes[x] = 0;
+					}
+					double r = random::U<double>();	// uniformly distributed random number betwee 0 and 1
+					if (r < pRes[x])
+					{
+						pRes[x] = 1;
+					}
+					else
+					{
+						pRes[x] = 0;
+					}
+				}
+			}
+			return res;
+		}
+
+		void CRBM::debug() {
+			std::cout << "Weight - rows: " << m_vpNeuronLayers[1]->getWeights().rows << " cols: " << m_vpNeuronLayers[1]->getWeights().cols << std::endl;
+
+			std::cout << "Positive H mean - rows: " << m_positiveHMean.rows << " cols: " << m_positiveHMean.cols << std::endl;
+			std::cout << "Positive H sample - rows: " << m_positiveHSample.rows << " cols: " << m_positiveHSample.cols << std::endl;
+			std::cout << "Negative H mean - rows: " << m_negativeHMean.rows << " cols: " << m_negativeHMean.cols << std::endl;
+			std::cout << "Negative H sample - rows: " << m_negativeHSample.rows << " cols: " << m_negativeHSample.cols << std::endl;
+			std::cout << "Negative V mean - rows: " << m_negativeVMean.rows << " cols: " << m_negativeVMean.cols << std::endl;
+			std::cout << "Negative V sample - rows: " << m_negativeVSample.rows << " cols: " << m_negativeVSample.cols << std::endl;    
+		}
+
+		void CRBM::sampleVisible(const Mat& values) {
+			m_negativeVMean = propagateDown(values);
+			m_negativeVSample = getBinomial(m_negativeVMean);
+		}
+
+		void CRBM::sampleHiddenPositive(const Mat& values) {
+			m_positiveHMean = propagateUp(values);		
+			m_positiveHSample = getBinomial(m_positiveHMean);
+
+			/*for (int y = 0; y < sample.rows; y++) {
+				float* pRess = sample.ptr<float>(y);
+				for (int x = 0; x < sample.cols; x++)
+					std::cout << pRess[x] << std::endl;
+			}*/
+		}
+
+		void CRBM::sampleHiddenNegative(const Mat& values) {
+			m_negativeHMean = propagateUp(values);
+			m_negativeHSample = getBinomial(m_negativeHMean);
+		}
+
+		Mat CRBM::propagateUp(const Mat& values) {
+			m_vpNeuronLayers[0]->setNetValues(values); //set the visible layer values
+
+			m_vpNeuronLayers[1]->dotProd(m_vpNeuronLayers[0]->getValues()); //sigmoid(sum(visible * weights)+bias)
+
+			return m_vpNeuronLayers.back()->getValues();
+		}
+
+		Mat CRBM::propagateDown(const Mat& values){
+			m_vpNeuronLayers[0]->dotProdVis(values, m_vpNeuronLayers[1]->getWeights()); 
+
+			return m_vpNeuronLayers[0]->getValues();
+		}
+
+		void CRBM::gibbsHVH(const Mat& hiddenSample) {
+			sampleVisible(hiddenSample);
+			sampleHiddenNegative(m_negativeVSample);
+		}
+		/* This implementation of RBM uses single step contrastive divergence algorithm, called CD-1  */
+		void CRBM::contrastiveDivergence(const Mat& values, float learningRate) {
+			//-------POSITIVE PHASE--------------------
+			/*In the positive phase, the input sample “v” from the visible layer is “clamped” to the input layer, 
+			and then is propagated to the hidden layer. The result of the hidden layer activation is h.    */
+			sampleHiddenPositive(values);
+
+			//------NEGATIVE PHASE---------------------
+			/*In the negative phase, “h” from the hidden layer is propagated back to the visible layer with the 
+			new v, say v’. This is then propagated back to the hidden layer with activation result “h”    */
+			gibbsHVH(m_positiveHMean);
+
+			std::vector<double> test = m_negativeHSample;
+			for (int i = 0; i < m_vpNeuronLayers[1]->getNumNeurons(); i++) {
+				//std::cout << i << std::endl;
+				for (int j = 0; j < m_vpNeuronLayers[0]->getNumNeurons(); j++)
+				{
+					m_vpNeuronLayers[1]->getWeights().at<float>(j, i) +=
+						learningRate * (m_positiveHMean.at<float>(i, 0) * values.at<float>(j, 0) - m_negativeHMean.at<float>(i, 0) * m_negativeVSample.at<float>(j, 0))/4000; // divide
+				}
+				m_vpNeuronLayers[1]->getBiases().at<float>(i, 0) += learningRate * (m_positiveHSample.at<float>(i, 0) - m_negativeHMean.at<float>(i, 0))/4000; //divide
+			}
+			for (int i = 0; i < m_vpNeuronLayers[0]->getNumNeurons(); i++)
+			{
+				//std::cout << i << std::endl;
+				m_vpNeuronLayers[0]->getBiases().at<float>(i, 0) += learningRate * (values.at<float>(i, 0) * m_negativeVSample.at<float>(i, 0))/4000; //divide
+			}
+		}
+
+		Mat CRBM::reconstruct(const Mat& values) {
+			Mat h, temp;
+
+			h = propagateUp(values);
+			temp = propagateDown(h);
+			return temp;
+		}
+	}
+}
\ No newline at end of file
diff --git a/modules/DNN/RBM.h b/modules/DNN/RBM.h
new file mode 100644
index 00000000..9c899107
--- /dev/null
+++ b/modules/DNN/RBM.h
@@ -0,0 +1,44 @@
+#pragma once
+
+#include "types.h"
+#include "NeuronLayer.h"
+
+namespace DirectGraphicalModels {
+	namespace dnn {
+		class CRBM {
+			public:
+				DllExport CRBM(const std::vector<ptr_nl_t>& vpLayers);
+				DllExport CRBM(const CRBM&) = delete;
+				DllExport ~CRBM(void) = default;
+
+				DllExport bool operator=(const CRBM&) = delete;
+
+				DllExport void debug(); //printing out rows and cols just to be sure they are correct
+
+				DllExport Mat getBinomial(const Mat& mean); //binomial ddistribution
+				DllExport Mat propagateUp(const Mat& values); //hidden layer = sigmoid(sum(visible*weights)+bias)
+				DllExport Mat propagateDown(const Mat& values); //visible layer = sigmoid(sum(hidden*weights)+bias)
+
+				DllExport void sampleVisible(const Mat& values); //sample viible for negative phase
+				DllExport void sampleHiddenPositive(const Mat& values); //sample hidden for positive phase
+				DllExport void sampleHiddenNegative(const Mat& values); //sample hidden for negative phase
+
+				DllExport void gibbsHVH(const Mat& hiddenSample); //gibbs sampling(sample visible and then sample hidden negative phase)
+				DllExport void contrastiveDivergence(const Mat& values, float learningRate);
+				DllExport Mat reconstruct(const Mat& values);
+				//DllExport void gibbsHVH()
+
+			private:
+				std::vector<ptr_nl_t>    m_vpNeuronLayers;
+
+				Mat						 m_positiveHMean;
+				Mat						 m_positiveHSample;
+		
+				Mat						 m_negativeHMean;
+				Mat						 m_negativeHSample;
+
+				Mat						 m_negativeVMean;
+				Mat						 m_negativeVSample;
+		};
+	}
+}
\ No newline at end of file