bossbpm.m

have some details about MAP

function [h, h_oneterm, Lest] = bossbpm(X, y, s2, g2, p)
%
% function [h, h_oneterm, Lest] = bossbpm(X, y, s2, g2, p)
%
% Function for computing the MMSE estimate for the linear regression
% problem y = Xh + e, when the model order is unknown and for
% computing an estimate of the model order L. These estimates are
% denoted BPM (Bayesian Parameter estimation Method) and BOSS
% (Bayesian Order Selection Strategy).
%  
% See the article "Empirical Bayes Linear Regression with Unknown
% Model Order" by Y. Sel閚, E. G. Larsson in ICASSP 2007, Honolulu,
% Hawaii, USA.
% 
% INPUT:
%  X    - The full regressor matrix of dimensions N by n.
%  y    - The observation vector of length N.
%  s2   - The variance of the white Gaussian noise e (scalar).
%  g2   - The length n vector of the variances of the parameter
%         vector h (for equal variance, g2 = gamma2 * ones(n,1)).
%         Note that all elements in g2 must be > 0.
%  p    - The length n vector of prior probabilities
%         for the considered model orders [1, ..., n].
%  
% OUTPUT:
%  h    - The MMSE estimate of the parameter vector obtained by BPM
%  h_oneterm - The ML/LS estimate of the model of order Lest
%  Lest - The BOSS model order estimate
%

% Copyright (C) 2007  Yngve Sel閚
%
% 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閚, 14 February, 2007. [e-mail], http://www.it.uu.se/katalog/ys/
%
% Yngve Sel閚
% Systems & Control, Dept of Information Technology
% Uppsala University
% P O Box 337, SE 751 05 UPPSALA, Sweden
%

thr = 1e-6; % threshold for numerical errors

if s2 == 0
  warning('You have specified s2 = 0, i.e. that the data are noisefree.')
end
if abs(sum(p)-1) > thr
  warning('The prior probabilities in p do not sum to unity.')
end
if length(find(p < 0)) > 0
  error('Elements in p cannot be negative.')
end
if (length(p) ~= length(g2)) | (length(p) ~= size(X,2))
  error('The length of p, g2 and the width of X need all be the same.')
end

y = y(:);   % columnize y
N = length(y);
n = size(X,2);

H = zeros(n, n);

XtX = X(:,1)'*X(:,1);
Zinv = 1/(1/g2(1) + 1/s2*XtX);

for k = 1:n
  if p(k) == 0
    logpHi_y(k,1) = -Inf;
  else
    logPHi = log(p(k));
    Qinv = eye(N)/s2 - 1/(s2*s2)*X(:,1:k)*Zinv*X(:, 1:k)';
    logdetQi = N*log(s2) + log(det(XtX*diag(g2(1:k))/s2 + eye(k)));
    H(1:k, k) = diag(g2(1:k)) * X(:,1:k)'*Qinv*y; % E(h|y)
    logPy_Hi = -.5 * logdetQi -.5 * y'*Qinv*y; % neglect
                                               % 2*pi-constant
    logpHi_y(k,1) = logPy_Hi + logPHi; % neglect p(y)
  end
  if k ~= n % compute quantities for next order
    xtX = X(:,(k+1))'*X(:,1:k);
    XtX = [XtX, xtX'; xtX, X(:,(k+1))'*X(:,(k+1))];
    Zinv = [Zinv, zeros(k,1); zeros(1,k), 0] + ...
              [-Zinv*xtX'/s2; 1]*[-xtX*Zinv/s2, 1] / ...
           (X(:,(k+1))'*X(:,(k+1))/s2 + 1/g2(k+1) - xtX*Zinv* ...
            (xtX)'/(s2*s2));
  end
end

% Find the MMSE-estimate
h = zeros(n,1);
h = H * exp(logpHi_y-max(logpHi_y)) / sum(exp(logpHi_y-max(logpHi_y)));

% Find the estimate with max P(Hi|y)
[junk,Lest] = max(logpHi_y);
h_oneterm = zeros(n,1);
h_oneterm = H(:, Lest);