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. ys@it.uu.se, 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 errorsif s2 == 0  warning('You have specified s2 = 0, i.e. that the data are noisefree.')endif abs(sum(p)-1) > thr  warning('The prior probabilities in p do not sum to unity.')endif length(find(p < 0)) > 0  error('Elements in p cannot be negative.')endif (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.')endy = y(:);   % columnize yN = 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));  endend% Find the MMSE-estimateh = 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);