I have been reading the book Pattern Recognition and Machine Learning (Bishop) for a while, and recently I came across this figure, which was created using Bernoulli mixture model on the MNIST dataset:
I figured it would be fun to code this, so I basically followed their algorithm:
Suppose you have N images of handwritten digits from 2 to 4. Let D be the number of pixels of the image (each MNIST image is 28 x 28, so D = 28 * 28)
You may want to model this image using a Bernoulli distribution for each pixel (remember a Bernoulli distribution with parameter 'mu' is just like flipping a damaged coin that has the probability
muof landing heads, and
1 - mufor landing tails). So you might toss a coin for each pixel, each with a different
mu, if it lands head you set the pixel to 1, otherwise you set it to 0. But with this model, you can easily see that each pixel is independent of each other and there is no way you can model handwritten digits with this.
Here is where the Bernoulli mixture model comes into play. Instead of using a Bernoulli for each pixel, we use a mixture of Bernoullis (that is, a weighted sum of Bernoullis), and this can be solved by using an algorithm called Expectation - Maximization. With this model, modelling digits suddenly becomes possible.
Of course this is just an informal treatment of the method, and for further readings you should definitely check out Chapter 9, Pattern recognition and machine learning, Section 9.2: Mixture of Bernoulli distribution
Let me summarize my results: I was able to recreate the figures in the book, and when I used 600 training images, I got 75% accuracy on the test set. But when I used all the training images that were from 2 to 4, the results were 90% (not bad for my first project, actually). Also the algorithm apparently worked very well with classifying the number 4, but it made many errors at number 2 (overfitting, anyone?).
However the code is not vectorized, and I fear that my coding style is not good. So any feedback would be great. Here is my code:
function [Labels] = ReadLabelsMNIST(filename) % Read the labels of the MNIST dataset % Written by Dang Manh Truong fp = fopen(filename, 'rb'); assert(fp ~= -1, ['Could not open', filename, '']); magic = fread(fp, 1, 'int32', 0, 'ieee-be'); assert(magic == 2049, ['Bad magic number in ', filename, '']); numLabels = fread(fp, 1, 'int32', 0, 'ieee-be'); Labels = fread(fp,inf, 'unsigned char'); assert(size(Labels,1) == numLabels, 'Mismatch in label count'); fclose(fp); end
Loading MNIST dataset:
function [images, Labels, numRows, numCols] = LoadMNIST(SelectedNumbers, type, numImages) % Preprocessing code for Bernoulli mixture model % Load a number of random images from the MNIST dataset % All the digits are within the SelectedNumbers row array % Written by Dang Manh Truong % type = 1 : Train set. type = 2 : Test set assert( (type == 1) | (type == 2), 'Type = 1 or Type = 2'); if type == 1 Labels = ReadLabelsMNIST('train-labels.idx1-ubyte'); filename = 'train-images.idx3-ubyte'; else Labels = ReadLabelsMNIST('t10k-labels.idx1-ubyte'); filename = 't10k-images.idx3-ubyte'; end fp = fopen(filename, 'rb'); assert(fp ~= -1, ['Could not open ', filename, '']); magic = fread(fp, 1, 'int32', 0, 'ieee-be'); assert(magic == 2051, ['Bad magic number in ', filename, '']); [~] = fread(fp, 1, 'int32', 0, 'ieee-be'); numRows = fread(fp, 1, 'int32', 0, 'ieee-be'); numCols = fread(fp, 1, 'int32', 0, 'ieee-be'); % Find 'numImages' random images that are from 2 to 4 % TODO: Replace 'find' with logical indexing % Index = find(2 <= Labels & Labels <= 4); Index = find(sum(bsxfun(@eq,Labels, SelectedNumbers),2) > 0); s = RandStream('mt19937ar','Seed',0); Permuted = randperm(s,size(Index,1)); % Permuted = randperm(size(Index,1)); Permuted = Permuted(1: numImages); Index = sort(Index(Permuted)); Labels = Labels(Index); images = zeros(numRows, numCols, numImages); prev = 0; ImageSize = numCols * numRows; for i = 1 : numImages % Ignore unneeded images fread(fp, (Index(i) - prev - 1) * ImageSize , 'unsigned char'); % Read image Temp = fread(fp,ImageSize, 'unsigned char'); Temp = reshape(Temp, numRows, numCols); images(:,:,i) = Temp; prev = Index(i); end fclose(fp); images = permute(images,[2 1 3]); % Reshape to #pixels x #examples images = reshape(images, size(images, 1) * size(images, 2), size(images, 3)); % Convert to double and rescale to [0,1], then binarize the images images = double(images) / 255; images(images < 0.5) = 0; images(images >= 0.5) = 1; end
Function to help find cluster for each image. Used in testing phase:
function [Cluster] = GetClusterBMM(images,mu, K) % Input : images : The input images (N x D) % : mu : Parameters for Bernoulli distribution for each pixel (K x D) % : phi: Mixing coefficients (K x 1) % : K: Number of mixtures % Output: Cluster (N x 1) The clusters most likely to be associated with % each image N = size(images,1); ClusterSum = zeros(N,K); Temp1 = mu; Temp2 = 1 - Temp1; for n = 1 : N for k = 1 : K % In: http://blog.manfredas.com/expectation-maximization-tutorial/ % They used the sum of mu's, but to be honest I don't know why % Anyway it only gives about 80% accuracy, while mine gives 90% %ClusterSum(n,k) = sum(Temp1(k,images(n,:) == 1)) + sum(Temp2(k,images(n,:) == 0)); ClusterSum(n,k) = prod(Temp1(k,images(n,:) == 1)) * prod(Temp2(k,images(n,:) == 0)); end end [~,Cluster] = max(ClusterSum,,2);
Training the model:
function [Correct, MisClassified] = TestBMM(X, TestX, mu, Labels, TestLabels) % Testing phase for Bernoulli mixture model % Written by Dang Manh Truong % The parameters here closely resemble those in the book Pattern % recognition and machine learning (Bishop), chapter 9 % N : Number of data points % N': Number of test points % K : Number of mixtures % D : Dimension of each data points % Input: X (N x D) Train data % : TestX (N' x D) Test data % : mu (K x D) Bernoulli parameters learned from training phase % : numTestImages Number of test data needed % : Labels (N x 1) Train labels % : TestLabels (N' x 1) Test labels % Output: Correct: The number of times the algorithm get it right % : MisClassified(10,10) : The misclassification matrix % MisClassified(i,j) : The number of times that the digit 'i' is % misclassified as digit 'j'. Of course the diagonal is zero K = size(mu,1); N = size(X,1); numTestImages = size(TestX,1); Correct = 0; MisClassified = zeros(10,10); digitsInTheSameCluster = zeros(10,numTestImages); TrainClusters = GetClusterBMM(X,mu,K); % N x 1 TestClusters = GetClusterBMM(TestX,mu,K); for i = 1 : numTestImages for n = 1 : N if TestClusters(i) == TrainClusters(n) digitsInTheSameCluster(Labels(n),i) = digitsInTheSameCluster(Labels(n),i) + 1; end end [~, AssignedLabel] = max(digitsInTheSameCluster(:,i)); if AssignedLabel == TestLabels(i) Correct = Correct + 1; else MisClassified(TestLabels(i),AssignedLabel) = MisClassified(TestLabels(i), AssignedLabel) + 1; end end
Testing the model:
function [mu, phi, Res, effNum] = TrainBMM(X, mu, phi, Res, effNum) % Training phase for Bernoulli mixture model using Expecation-Maximization % Written by Dang Manh Truong % The parameters here closely resemble those in the book Pattern % recognition and machine learning (Bishop), chapter 9 % D: The dimension of each data points % N: The number of data points % K: The number of Bernoulli mixtures % Input: X (N x D) Data points to be processed (each row - a data point) % : mu(K x D) Bernoulli parameters for each mixture % : phi(K x 1) Mixing coefficients for each mixture % : Res(N x K) Responsibilities of each component (1 - K) given a data % point (1 - n) % : effNum(K x 1) Effective number of observations for each mixture % Output: The new values of mu, phi, Res and effNum % Most of the time only mu will be used % Size of each image. I don't want to pass these to the function % because their only purpose is to show the images numRows = 28; numCols = 28; N = size(X,1); K = size(phi,1); iterNum = 0; uniform = 1 / K; fprintf('E-M algorithm in progress. This may take a while.....\n'); while 1 % E-step % Equivalent unvectorized code: % for n = 1 : N % for k = 1 : K % Res(n,k) = 1; % for i = 1 : D % D = size(X,2) % if X(n,i) == 1 % Res(n,k) = Res(n,k) * mu(k,i); % else % Res(n,k) = Res(n,k) * (1 - mu(k,i)); % end % end % end % end % TODO: Vectorize this part ASAP!!!! for n = 1 : N for k = 1 : K Temp1 = mu(k,:); Temp2 = 1 - mu(k,:); Res(n,k) = prod(Temp1(X(n,:) == 1)) * prod(Temp2(X(n,:) == 0)); end end Res = bsxfun(@times, Res, phi'); % Divide by the denominator Sum = sum(Res,2); Sum(Sum == 0) = uniform; Res = bsxfun(@rdivide, Res, Sum); % M-step effNum = sum(Res,1); mu = Res' * X; mu = bsxfun(@rdivide, mu, effNum'); % Check for convergence iterNum = iterNum + 1; for k = 1 : K subplot(1,K,k) Result = reshape(mu(k,:), numRows, numCols); subimage(Result) end hold on pause(1) fprintf('Iteration %d \n',iterNum); if iterNum >= 10 break; end end fprintf('Press any key to continue\n\n\n'); pause close(gcf) end
The main function:
function  = BMM() % Bernoulli mixture model for classification of MNIST dataset % Based on Figure 9.10 in the book Pattern recognition and machine learning % Partly inspired by: http://blog.manfredas.com/expectation-maximization-tutorial/ % Written by Dang Manh Truong % I was able to reproduce the 3 pictures of digits 2,3 and 4, and the % results for classfication on the MNIST test set (2 to 4) were satisfying % with about 90% correct, not too bad for my first project. % But it's not a good score compared with state-of-the-art methods pause(1); fprintf('Bernoulli mixture model using Expectation - Maximization\n'); fprintf('to recreate Figure 9.10, chapter 9, Pattern recognition and machine learning\n'); % Change these lines if you wish SelectedNumbers = [2 3 4 ]; % The numbers that we care about in the dataset numTrainImages = 600; % For the numbers from 2-4: <= 17391 numTestImages = 3024; % For those from 2-4: <= 3024 rng(0,'twister'); % Step 1: Initialization % images: #pixels * #examples [images, Labels, numRows, numCols] = LoadMNIST(SelectedNumbers, 1, numTrainImages); N = numTrainImages; % N : The number of train images K = size(SelectedNumbers,2); % K : The number of mixtures D = numRows * numCols; % Dimension of each image phi = ones(K,1) * 1/K; % Mixing coefficients mu = (0.75-0.25) * rand(K,D) + 0.25 ; % Means of each components mu = mu ./ repmat(sum(mu,2),1,D); Res = zeros(N,K); % Res(k,n): Responsibility of component 'k' given data point X(n,:) effNum = zeros(K,1); % Effective number of data points associated with each component X = images'; % (N x D) Each row is an image % Step 2: Expectation - Maximization [mu, ~, ~, ~] = TrainBMM(X, mu, phi, Res, effNum); % Step 3: Testing [TestImages, TestLabels ,~, ~] = LoadMNIST(SelectedNumbers, 2,numTestImages); TestX = TestImages'; % Each image is in one row [Correct, MisClassified] = TestBMM(X, TestX, mu, Labels, TestLabels); fprintf('Correct: %f percents \n',100 * Correct / numTestImages); fprintf('The misclassification matrix: \n'); MisClassfied = MisClassified(SelectedNumbers, SelectedNumbers) % The MisClassified matrix when I used all 17931 train images (that are % from 2 to 4) to test all 3024 test images (again, from 2 to 4): % 0 130 48 % 84 0 21 % 17 1 0 % The algorithm appears to correctly labels all the digits 4, but fails % spectacularly at the digit 2.
To run this program you need the MNIST dataset, available here: http://yann.lecun.com/exdb/mnist/ .Just type "BMM()" in the Matlab command line, and the program will run.
Here are the 3 resulting clusters: