% Matlab script to reproduce Figure 1 of 
% Gribonval, R., "Should penalized least squares regression be interpreted
% as Maximum A Posteriori estimation", Signal Processing, IEEE Transactions
% on, 2011, to appear.


%% Set the parameters of a Bernoulli-Gaussian distribution, see Eq. III.1
p=0.99;      % probability of zero coefficient
sigma2_1=10;  % variance of nonzero coefficients
a=((1-p)/p)*sqrt(1/(sigma2_1+1));
b=1-1/(sigma2_1+1);
c=sigma2_1/(sigma2_1+1);

%% Figure 1 (b): plot prior and evidence, using analytic expression

% The variable y
dy = 0.01; 
y=0:dy:15; 
ysym = [-y(end:-1:2) y];
% The evidence p_Y(y)
pY   = p*exp(-y.^2/(2*(0+1)))/sqrt(2*pi*(0+1))+(1-p)*exp(-y.^2/(2*(sigma2_1+1)))/sqrt(2*pi*(sigma2_1+1));
pYsym= [pY(end:-1:2) pY];
% The prior p_X(y)
pX   = (1-p)*exp(-y.^2/(2*sigma2_1))/sqrt(2*pi*sigma2_1);
pXsym= [pX(end:-1:2) 0.5*pY(1) pX(2:end)];
% Plot on Figure 1 b
h=figure(1);
set(h,'Position',[15         326        1000         400]);
subplot(1,2,2);
hold off;
plot(ysym,-log(pXsym),'b:');
hold on;
plot(ysym,-log(pYsym),'r-.');
%legend('\phi_{MAP}(x) = -log p_X(x)','-log p_Y(x)','Location','North');
xlabel('x');
title(['(b) Penalty functions (p=' num2str(p) ' \sigma_1^2=' num2str(sigma2_1) ')']);

%% Build the MMSE estimator function as in Eq. (41), p.11 of reference
%% [12], draft version:
% J.S. Turek, I. Yavneh, M. Protter, M. Elad, "On MMSE and MAP Denoising
% Under Sparse Representation Modeling Over a Unitary Dictionary",
% submitted to Applied and Computational Harmonic Analysis. Technion
% University, 2010.

% The MMSE estimator
psiM = c*y./(1+exp(-b*y.^2)/a);
psiMsym = c*ysym./(1+exp(-b*ysym.^2)/a);
% Plot on Figure 1 a
subplot(1,2,1);
plot(ysym,psiMsym,'b');
xlabel('y');
title(['(a) MMSE estimator and its inverse (p=' num2str(p) ' \sigma_1^2=' num2str(sigma2_1) ')']);

%% Build the associated penalty function
% From Eq. (II.2), compute the derivative of  
% dlogq(y) := \nabla log p_Y(y) = \psi_MMSE(y)-y
dlogqsym = psiMsym-ysym;  
% From the definition logq(y) := log p_Y(y)
logqsym  = log(pYsym);
% From Eq. (II.7)  \phi[\psi_MMSE(y)] = -||y-\psi_MMSE(y)||_2^2-\log p_Y(y)
phipsisym= -(dlogqsym.^2)/2 -logqsym;
% Plot on Figure 1 b
subplot(1,2,2);
hold on;
plot(psiMsym,phipsisym,'b');
hold off;
h = legend('\phi_{MAP}(x) = -log p_X(x)','-log p_Y(x)','\phi_{MMSE}(x)','Location','North');
%set(h,'Interpreter','latex');
%set(h,'String',{'\varphi_{MAP}(x) = -log p_X(x)','-log p_Y(x)','\varphi_{MMSE}(x)'});

%% Build the inverse function
dx=dy;
x=0:dx:max(psiM);
for i=1:length(x)
     idx(i) = min(find(psiM>=x(i)));
     invpsiM(i)= y(idx(i));
end

% Plot on Figure 1a
subplot(1,2,1);
hold on;
plot([-x(end:-1:2) x],[-invpsiM(end:-1:2) invpsiM],'r-.');
hold off;
h=legend('\psi_{MMSE}(y)','\psi_{MMSE}^{-1}(y)','Location','North');