empbayes_pen.m

have some details about MAP

function gammaest = empBayes_pen(hLS, lam, Lmat)
% function gammaest = empBayes_pen(hLS, lam, Lmat)
%  
% Function for estimating the coefficient variances gamma for a linear
% regression equation. This function uses the "Estimation by
% Penalization"-approach described in Section 4.2 of "Empirical
% Bayes Linear Regression with Unknown Model Order" by Y. Sel閚,
% E. G. Larsson (submitted to ICASSP 2007). 
%
% INPUT:
%  hLS  - The full model LS/ML coefficient estimate.
%  lam  - The smoothing parameter lamda (see eq (18) in the above reference).
%  Lmat - The difference matrix (default, second order).
%
% 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

n = length(hLS); % full order

if nargin < 3 % if no Lmat is given, then construct second order
              % difference matrix
  Lmat = diag(ones(n,1)) - 2*diag(ones(n-1,1),1) + diag(ones(n-2,1), 2);
  Lmat = Lmat(1:n-2,:);
end

% Start with closed from version, which may give negative variance
% estimates. If it does, then proceed with constrained quadratic
% programming. 
gammaest = inv(eye(n) + lam*Lmat.'*Lmat) * abs(hLS).^2;
if length(find(gammaest < 0)) % if negative elements in gammaest
  gammaest = quadprog(eye(n) + lam*Lmat.'*Lmat, -abs(hLS).^2, ...
                      -eye(n), zeros(n,1), [], [], [], [], abs(hLS).^2, ...
                      opt);
end