function [P_RS, P_S] = light_estimation(EVEN_S_DIST)
% Demo MAP, ML, Bayes methods in estimation
% WW/BME/WU     02/08/2004
rand('state',0);
S = 0:20;
R = 0:20;
Reading = [3 10];

d = (max(S)-min(S))/8;
d2 = d/4;
%Generate P(R|S)
gaussian = 1/(sqrt(2*pi)*d)*exp(-((S-mean(S)).^2)/(2*d^2));
gaussian = gaussian/sum(gaussian)/2;
two_humps = 1/(sqrt(2*pi)*d)*exp(-((S-mean(S)/2).^2)/(2*d^2)) + 1/(sqrt(2*pi)*d)*exp(-((S-mean(S)*1.5).^2)/(2*d^2))/2;
two_humps = two_humps/sum(two_humps)/2;
rest = 1 - gaussian - two_humps;
rest = repmat(rest, length(R)-2, 1);
random = rand(length(R)-2,length(S));
rand_sum = sum(random);
rand_sum = repmat(rand_sum, length(R)-2, 1);
P_RS([1:Reading(1)-1 Reading(1)+1:Reading(2)-1 Reading(2)+1:length(R)],:)  = random./rand_sum.*rest;
P_RS(Reading(1),:) = gaussian;
P_RS(Reading(2),:) = two_humps;
if EVEN_S_DIST
    P_S = 1/length(S)*ones(1, length(S));
else
    P_S = 1:length(S);
    P_S = P_S/sum(P_S);
end
P_R = P_RS * P_S';
%Compute P(S|R)
for i = 1:length(R)
    for j = 1:length(S)
        P_SR(i,j) = P_RS(i,j)*P_S(j)/P_R(i);
    end
end


for r = Reading
    fprintf('Reading = %d:\n',r);
    figure
    hold on
    plot(S, P_RS(r,:), 'LineWidth',2)
    plot(S, P_SR(r,:),'r','LineWidth',2)
    legend('P(R|S)  ML', 'P(S|R)  MAP')
    title(sprintf('Reading = %d', r))

    %Compute Bayes error
    for i = 1:length(S)
        Bayes1(i) = sum((S - S(i)).^2.*P_SR(r,:));
        Bayes2(i) = sum(abs(S-S(i)).*P_SR(r,:));
    end

    %Normalize the error for easy comparision
    Bayes1 = (Bayes1-min(Bayes1))/(max(Bayes1)-min(Bayes1));
    Bayes2 = (Bayes2-min(Bayes2))/(max(Bayes2)-min(Bayes2));

    figure
    hold on
    plot(S, Bayes1, 'LineWidth',2)
    plot(S, Bayes2, 'r','LineWidth',2)
    legend('Bayes1 (S-S_b_a_y_e_s)^2', 'Bayes2 |S-S_b_a_y_e_s|')
    title(sprintf('Reading = %d', r))
    
    %Print answer
    [m,MAP] = max(P_SR(r,:));
    [m,ML] = max(P_RS(r,:));
    [m,B1] = min(Bayes1);
    [m,B2] = min(Bayes2);
    mean_p_sr = sum(P_SR(r,:).*S);
    for i = 2:length(S)-1
        P_pre = sum(P_SR(r,1:i-1));
        P_post= sum(P_SR(r,i+1:end));
        P_diff(i-1) = abs(P_post - P_pre);
        [m,median_p_sr] = min(P_diff);
    end
    fprintf('MAP    :%d\n', MAP-1);
    fprintf('ML     :%d\n', ML-1);
    fprintf('Bayes1 :%d    mean    :%d\n', B1-1, round(mean_p_sr));
    fprintf('Bayes2 :%d    median  :%d\n', B2-1, median_p_sr);
end
