diff --git a/source/Rbm.cpp b/source/Rbm.cpp
index 0faa76b..e284ba3 100644
--- a/source/Rbm.cpp
+++ b/source/Rbm.cpp
@@ -252,23 +252,20 @@ void Rbm::cd_hinton(arma::mat const &v_data, arma::mat &dw, arma::mat &dbh, arma
 void Rbm::cd_jens(arma::mat const &v_states, arma::mat &dw, arma::mat &dbh, arma::mat &dbv)
 {
 	arma::mat v_probs(v_states);
-	arma::mat h_states = toHiddenProbs(v_states);
-	arma::mat h_probs;
-	
+	arma::mat h_probs = toHiddenProbs(v_states);
+	arma::mat h_states;
+
 	// Sample hidden
 	if (m_params.doGaussianHidden)
 	{
-		h_probs = h_states;
 		h_states = h_probs + arma::randn(h_probs.n_rows, h_probs.n_cols);
 	}
 	else if (m_params.doRaoBlackwell)
 	{
-		h_probs = toHiddenProbs(v_states);
 		h_states = h_probs;
 	}
 	else
 	{
-		h_probs = toHiddenProbs(v_states);
 		h_states = sample(h_probs);
 	}
 	
@@ -313,8 +310,7 @@ void Rbm::train(arma::mat const &batch, IListener* pListener)
 	double dProgress = 100.0/(batch.n_rows*m_params.numEpochs);
 	double progress = 0;
 	int lastProgress = -100;
-	int batchRowIndex = 0;
-	
+
 	arma::mat dbv(arma::zeros(1, m_bv.n_cols));
 	arma::mat dbh(arma::zeros(1, m_bh.n_cols));
 	arma::mat dwhv(arma::zeros(m_whv.n_rows, m_whv.n_cols));
@@ -322,40 +318,46 @@ void Rbm::train(arma::mat const &batch, IListener* pListener)
 	arma::mat inc_bv(arma::zeros(1, m_bv.n_cols));
 	arma::mat inc_bh(arma::zeros(1, m_bh.n_cols));
 
-	int trainingSizeRemain = batch.n_rows;
-
 	bool shouldAbort = false;
-	while (trainingSizeRemain && !shouldAbort)
+	// One epoch = one full pass over the (freshly reshuffled) data, mini-batch
+	// by mini-batch. Epochs must be the outer loop: nesting mini-batches on
+	// the outside would run all numEpochs steps on one fixed chunk before ever
+	// moving to the next, biasing training toward whatever data comes last.
+	for (int epoch=0; epoch < m_params.numEpochs && !shouldAbort; epoch++)
 	{
-		int miniBatchSizeActual = std::min(m_params.miniBatchSize, trainingSizeRemain);
-		arma::mat miniBatch = batch.rows(batchRowIndex, batchRowIndex+miniBatchSizeActual-1);
-		trainingSizeRemain -= miniBatchSizeActual;
-		batchRowIndex += miniBatchSizeActual;
-		int numcases = std::min(m_params.miniBatchSize, (int)batch.n_rows);
+		arma::mat shuffledBatch = arma::shuffle(batch);
+		int trainingSizeRemain = shuffledBatch.n_rows;
+		int batchRowIndex = 0;
 
-		arma::mat v_states(miniBatch);
+		while (trainingSizeRemain && !shouldAbort)
+		{
+			int miniBatchSizeActual = std::min(m_params.miniBatchSize, trainingSizeRemain);
+			arma::mat miniBatch = shuffledBatch.rows(batchRowIndex, batchRowIndex+miniBatchSizeActual-1);
+			trainingSizeRemain -= miniBatchSizeActual;
+			batchRowIndex += miniBatchSizeActual;
+			int numcases = std::min(m_params.miniBatchSize, (int)batch.n_rows);
+
+			arma::mat v_states(miniBatch);
+
+			// Create hidden layer base on training data
+			if (m_params.doSampleBatch)
+			{
+				// When the hidden units are being driven by data, always use stochastic binary states
+				v_states = sample(miniBatch);
+			}
 
-		// Create hidden layer base on training data
-		if (m_params.doSampleBatch)
-		{
-			// When the hidden units are being driven by data, always use stochastic binary states
-			v_states = sample(miniBatch);
-		}
-		
-		for (int epoch=0; epoch < m_params.numEpochs; epoch++)
-		{
 			// Contrastive divergence learning: calculate gradients
 			cd(v_states, dwhv, dbh, dbv);
-			
+
 			// Adjust weight and biases
 			inc_bv = m_params.momentum*inc_bv + m_params.learningRate/numcases*dbv;
 			inc_bh = m_params.momentum*inc_bh + m_params.learningRate/numcases*dbh;
 			inc_whv = m_params.momentum*inc_whv + m_params.learningRate*(dwhv/numcases - m_params.weightDecay*m_whv);
-			
+
 			m_bv += inc_bv;
 			m_bh += inc_bh;
 			m_whv += inc_whv;
-			
+
 			progress += dProgress*miniBatchSizeActual;
 			status.progress = (int)(progress + 0.5);
 
@@ -363,7 +365,7 @@ void Rbm::train(arma::mat const &batch, IListener* pListener)
 			if (status.progress != lastProgress)
 			{
 				lastProgress = status.progress;
-				
+
 				// Calculate error
 				status.err = rms_error_accu(miniBatch - toVisibleProbs(toHiddenProbs(v_states)));
 				if (pListener)
@@ -375,10 +377,10 @@ void Rbm::train(arma::mat const &batch, IListener* pListener)
 					}
 				}
 			}
-			
-		} // Number of epochs
-			
-	} // number of mini batches
+
+		} // mini batches
+
+	} // epochs
 
 	// Update final status
 	status.err_total = rms_error_accu(batch - toVisibleProbs(toHiddenProbs(batch)));