function [sample, probs] = rec (obs, digits, noise)
% Prepare the matrix p_n(X,Y) which is the probability that a pixel
% with value X in the digit template becomes a pixel with value Y after
% adding noise.
% convolution is an easy way to get copies of noise in shifted in a matrix
p_n = conv2 (eye(4), noise);
% now we threshold values into the 1..4 range
p_n = [sum(p_n(:,1:4),2), p_n(:,5:6), sum(p_n(:,7:10),2)];
% we use log-probabilities for more stability (and more speed too)
logp_n = log (p_n);
n_rows = size(digits, 1);
n_cols = size(digits, 2);
n_digits = size(digits, 3);
n_tests = size(obs, 3);
% output variables, for each test case:
% the most probable digit
sample = zeros (n_tests, 1);
% and probability per digit
probs = zeros (n_tests, n_digits);
% for each test case
for c = 1:n_tests
logp_d = zeros (1, n_digits);
% and for each digit
for d = 1:n_digits
% logp = 0 means a probability of 1
logp = 0;
% and in each row
for i = 1:n_rows
% for each pixel of it
for j = 1:n_cols
% the probability that the this digit evolves into the
% given observation, is the product of all per-pixel
% probabilities. for each pixel, the total probability
% that the pixel turned from X to Y after adding noise,
% is indeed p_n(X,Y). since we are using log-probabilities,
% we add, instead of multiply.
logp = logp + logp_n(digits(i,j,d), obs(i,j,c));
end
end
% fill in
logp_d(d) = logp;
end
% ok, now we computed all the likelihoods, i.e., the probability that
% digit d generated the observation, for each d. now to assign
% responsibilities to each digit, we need to normalize the probabilities,
% to add to 1. one of the digits should be responsible afterall.
% to do the normalization, one way is to simply convert from
% log-probabilities to probabilities, and divide by the sum of them (to
% make them sum to 1). we take another step before that, scale up the
% log-probabilities such that one of them becomes zero. this makes the
% following computation more stable, as it doesn't pressure the exponent.
logp_d = logp_d - max (logp_d);
p_d = exp(logp_d);
p_d = p_d / sum(p_d);
% now p_d is the desired probability distribution
% save it
probs(c,:) = p_d;
% moreover, main suspect will be the one with the highest
% probability/responsibility.
[best, sample(c)] = max (p_d);
end