empbayes_par.m

have some details about MAP

function gammaest = empBayes_par(Psi, Psiq, hLS)
% function gammaest = empBayes_par(Psi, Psiq, hLS)
%  
% Function for estimating the coefficient variances gamma for a linear
% regression equation. This function uses the "Estimation by
% Parameterization"-approach described in Section 4.1 of "Empirical
% Bayes Linear Regression with Unknown Model Order" by Y. Sel閚,
% E. G. Larsson, ICASSP 2007. 
%
% INPUT:
%  Psi - The matrix Psi in eq (15) of the above reference (we have
%  parameterized gamma such that gamma = Psi*alpha).
%  Psiq - The matrix diag[q_1, ..., q_n] * Psi in eq (17) of the
%         above reference.
%  hLS  - The full model LS/ML coefficient estimate.
%  
% OUTPUT:
%  gammaest - The estimated variances for the coefficient vector h.

% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License
% as published by the Free Software Foundation; either version 2
% of the License, or (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
% MA  02110-1301, USA.
%

% Yngve Sel閚, 3 October 2006, [e-mail], www.it.uu.se/katalog/ys/
  
opt=optimset('Largescale','off'); % for the linear program below

% Standard LS estimate of the alpha vector (see eq (15)). 
alpha0 = Psiq\[abs(hLS).^2];

% With the above reparameterization, gammaest = Psi * alpha0. If this
% gives negative gamma estimates, linear programming with
% constrained gamma (always > 0) must be used:
if length(find(Psi*alpha0<0)) > 0
  % Quadratic minimization with constraint rho = Psi*alpha > 0
  n = length(hLS);
  alpha = quadprog(Psiq'*Psiq, -Psiq'*[abs(hLS).^2], -Psi, ...
                   zeros(n,1), [], [], [], [], alpha0, opt);
  gammaest = Psi*alpha;
else
  gammaest = Psi*alpha0;
end