1 % Configuration
   2 % -------------
   3 nv = 784;   % number of visible units
   4 nh = 400;   % number of hidden units
   5 mb = 100;   % size of the minibatch
   6 it = 100;   % number of iterations (WARNING: set to 10000 or higher to produce the final results)
   7 
   8 % Initialize data
   9 % ---------------
  10 data = load('mnist.mat');       % Load the dataset
  11 X = single(data.x)/255.0 > 0.5; % Preprocess the dataset and convert it into a matrix of binary values
  12 W = normrnd(0,0.001,[nv,nh]);   % Weight matrix
  13 a = zeros(1,nv);                % Bias vector for visible units
  14 b = zeros(1,nh);                % Bias vector for hidden units
  15 
  16 % Algorithm
  17 % ---------
  18 for i=1:it,
  19 	
  20 	% Take a minibatch from the dataset
  21 	P = ceil(unifrnd(0,10000,[mb,1]));
  22 	v = X(P(:),:);
  23 	
  24 	% Perform three steps of Gibbs sampling
  25 	h = % TO COMPLETE
  26 	V = % TO COMPLETE
  27 	H = % TO COMPLETE
  28 	
  29 	% Compute the directions of parameter updates
  30 	dW = % TO COMPLETE
  31 	da = % TO COMPLETE
  32 	db = % TO COMPLETE
  33 	
  34 	% Update the parameters
  35 	W = W + 0.001 * dW;
  36 	b = b + 0.001 * db;
  37 	a = a + 0.001 * da;
  38 	
  39 	% Print iterations and some information on the parameters of the RBM and the reconstruction error
  40 	sprintf('i=%5d   std(W)=%.3f   mean(a)=%.3f   mean(b)=%.3f   e=%.3f\n',i,std(W(:)),mean(a(:)),mean(b(:)),mean(mean(abs(v-V))))
  41 end
  42 
  43 % Display the weights and the reconstructed digits
  44 % ------------------------------------------------
  45 wdisp = reshape(permute(reshape(W,[28,28,20,20]),[2,4,1,3]),[28*20,28*20]); % Create a human-readable view of matrix W
  46 wdisp = wdisp - min(wdisp(:));
  47 wdisp = wdisp / max(wdisp(:));
  48 imwrite(wdisp,'w.png')
  49 
  50 xdisp = reshape(permute(reshape(V,[10,10,28,28]),[4,2,3,1]),[28*10,28*10]);
  51 xdisp = xdisp - min(xdisp(:));
  52 xdisp = xdisp / max(xdisp(:));
  53 imwrite(xdisp,'x.png')