simulation.m

have some details about MAP

% This is a simulation script for testing the methods used in the
% article "Empirical Bayes Linear Regression with Unknown Model
% Order" by Y. Sel閚, E. G. Larsson. 

% 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/

clear all 

randomstate = sum(100*clock);
randn('state', randomstate);
rand('state', randomstate);

% simulation parameters

n = 30   % max nr of peaks to consider
Nt = 50; % length of data sequence
Nmonte = 10000; % number of Monte Carlo simulations

correl = 0;  % amount of correlation between regressors: 
           % X(:,k) = (correl*z + randn(Nt,1)) / sqrt(correl^2+1)

s2_db = 10:-2.5:-20;  % standard deviation of noise in dB (relative
                      % to 1)
s2 = 10.^(s2_db/10);  % noise standard deviation

% Statistics:
% MSEs:
MSE_BIC = zeros(length(s2), 1);
MSE_multiBIC = zeros(length(s2), 1);
MSE_AIC = zeros(length(s2), 1);
MSE_AICc = zeros(length(s2), 1);
MSE_Bayes = zeros(length(s2), 1);
MSE_Bayes_oneterm = zeros(length(s2), 1);
Bayes_beatsIC = zeros(length(s2), 1);
MSE_empBayesSmooth = zeros(length(s2), 1);
MSE_empBayesSmooth_oneterm = zeros(length(s2), 1);
empBayesSmooth_beatsIC = zeros(length(s2), 1);
MSE_empBayesSmoothF = zeros(length(s2), 1);
MSE_empBayesSmoothF_oneterm = zeros(length(s2), 1);
empBayesSmoothF_beatsIC = zeros(length(s2), 1);
MSE_empBayesSmoothD = zeros(length(s2), 1);
MSE_empBayesSmoothD_oneterm = zeros(length(s2), 1);
empBayesSmoothD_beatsIC = zeros(length(s2), 1);
MSE_empBayesSmoothGCV = zeros(length(s2), 1);
MSE_empBayesSmoothGCV_oneterm = zeros(length(s2), 1);
empBayesSmoothGCV_beatsIC = zeros(length(s2), 1);

% correct order
corrBIC = zeros(length(s2), 1);
corrAIC = zeros(length(s2), 1);
corrAICc = zeros(length(s2), 1);
corrBayes = zeros(length(s2), 1);
corrempBayesSmooth = zeros(length(s2), 1);
corrempBayesSmoothF = zeros(length(s2), 1);
corrempBayesSmoothD = zeros(length(s2), 1);
corrempBayesSmoothGCV = zeros(length(s2), 1);

% histograms of selected orders
ordBIC = zeros(length(s2),n);
ordAIC = zeros(length(s2),n);
ordAICc = zeros(length(s2),n);
ordBayes = zeros(length(s2),n);
ordempBayesSmooth = zeros(length(s2),n);
ordempBayesSmoothF = zeros(length(s2),n);
ordempBayesSmoothD = zeros(length(s2),n);
ordempBayesSmoothGCV = zeros(length(s2),n);

% For Empirical Bayes by Parameterization (Section 4.1). 
% True variance profiles:
v = exp(-0.4*(1:n)'); % exponential decay
v = v/sqrt(v'*v) * sqrt(n); % norm == n
v2 = exp(0.4*(1:n)'); % exponential increase
v2 = v2/sqrt(v2'*v2) * sqrt(n); % norm == n
Psi = [ones(n,1), v, v2];
r = size(Psi,2);
Psiq = Psi .* repmat((n:-1:1)'/n, 1, r);

% For Empirical Bayes by Parameterization (Section 4.1). 
% Fourier-based variance profiles:
kF = [1,1]; % [number of cos-terms, number of sin-terms] (normally "==")
PsiF = ones(n,1); % constant term
for k = 1:kF(1)
  PsiF = [PsiF, cos(k*2*pi/n*(0:n-1)')];
end
for k = 1:kF(2)
  PsiF = [PsiF, sin(k*2*pi/n*(0:n-1)')];
end
PsiFq = PsiF .* repmat((n:-1:1)'/n, 1, sum(kF)+1);

% Second order differentiation matrix for smoothing
Lmat = diag(ones(n,1)) - 2*diag(ones(n-1,1),1) + diag(ones(n-2,1), 2);
Lmat = Lmat(1:n-2,:);
Lo = 2; % order of difference equation

% GSVD for GCV (see [Hansen])
[U_,V_,X_,S_,M_]= gsvd(eye(n),Lmat);
gsv = diag(S_(1:n-Lo,1:n-Lo)) ./ diag(M_);
d0mat = eye(n) - U_*U_';
lamvec2 = [0.1, (0.5:.5:100),110:10:200, 300:100:1000, 5000, 10000]; 
nolam2 = length(lamvec2);

opt=optimset('Largescale','off');

tic
for nm = 1:Nmonte
  if (rem(nm,10) == 0) % print out iteration number
    nm
  end
  for nsnr = 1:length(s2) % for each noise level
      z = randn(Nt,1); % correlation vector
      % construct regressor matrix:
      X = (randn(Nt,n) + correl*repmat(z,1,n)) / sqrt(1+correl^2);
      Hest = zeros(n,n); % store coefficient estimates of different
                         % orders
      L = ceil(rand * n); % true number of coefficients, between 1 and n
      G = diag(Psi * randn(size(Psi,2),1).^2); % covariance
                                               % matrix of h
      G = G / sqrt(diag(G)' * diag(G));           % normalize G
      h = sqrt(G)*[randn(L,1);zeros(n-L,1)]; % true coeff vector
      e = sqrt(s2(nsnr))*randn(Nt,1); % generate Gaussian white noise
      y = X*h + e; % noisy signal

      % Find LS/ML solutions for AIC and BIC:
      for k = 1:n % for each order
         Hest(1:k, k) = X(:,1:k) \ y; % estimate filter coeff
         s2hat(k) = sum(abs(y - X*Hest(:, k)).^2) / Nt; % noisevar
                                                        % ML estimate
      end

      % Decide lower threshold for variance profile:
      bFac = 0.1;
      b = max(Hest(:,n).^2)*bFac;

      % unbiased estimate of s2hat (for largest order n)
      s2hat_ub = s2hat(n)*Nt/(Nt-n);
      
      % Moment-based peak variance estimation using estimated
      % variance profile combination from true set
      rho = empBayes_par(Psi, Psiq, Hest(:,n));
      
      % Implement the threshold that variances >= b:
      rho(find(rho<b)) = b; 
      
      % Moment-based peak variance estimation using estimated
      % variance profile from Fourier based Psi-matrix

      rhoF = empBayes_par(PsiF, PsiFq, Hest(:,n));
      
      % Implement the threshold that variances >= b:
      rhoF(find(rhoF<b)) = b; 

      % Tikhonov regularization type penalty for 2nd order
      % difference equation (min || res ||^2 + lam || diff(rho) ||^2
      
      lam = 3; % smoothness parameter
      
      rhoD = empBayes_pen(Hest(:,n), lam, Lmat);
      
      % Implement the threshold that variances >= b:
      rhoD(find(rhoD<b)) = b; 
                  
      % Find the appropriate smoothnes parameter lambda using GCV
      % from [Hansen] 1992
      Glam = sum([repmat(lamvec2,n-Lo,1)./(repmat(gsv.^2,1,nolam2) + ...
                                  repmat(lamvec2,n-Lo,1)) .* ...
                                  repmat(U_(:,1:n-Lo)'*abs(Hest(:, ...
                                                    n)).^2, 1, nolam2)].^2) ./ ...
          [sum(repmat(lamvec2,n-Lo,1)./(repmat(gsv.^2,1,nolam2) + ...
                                  repmat(lamvec2,n-Lo,1)))].^2;
      [slask,minind] = min(Glam);
      lamest = lamvec2(minind);
      
      rhoGCV = empBayes_pen(Hest(:,n), lamest, Lmat);
            
      % Implement the threshold that variances >= b:
      rhoGCV(find(rhoGCV<b)) = b; 
      

      % NOW, WE START ESTIMATING THE ORDERS AND THE COEFFICIENT
      % VECTORS:
      
      % The information criteria:
      BIC = Nt*log(s2hat) + log(Nt)*(1:n);
      AIC = Nt*log(s2hat) + 2*(1:n);
      AICc = Nt*log(s2hat) + 2*Nt ./ (Nt-2-(1:n)) .* ((1:n)+1);
      
      hmultiBIC = Hest * [exp(-.5 * (BIC-min(BIC))) / ...
          sum(exp(-.5 * (BIC-min(BIC))))].';
      
      % The selected orders:
      [slask, BICo] = min(BIC);
      [slask, AICo] = min(AIC);
      [slask, AICco] = min(AICc);
      
      % XICo contains [selected order]
      
      ordBIC(nsnr, BICo) = ordBIC(nsnr, BICo) + 1;
      ordAIC(nsnr, AICo) = ordAIC(nsnr, AICo) + 1;
      ordAICc(nsnr, AICco) = ordAICc(nsnr, AICco) + 1;
      
      MSE_BIC(nsnr) = MSE_BIC(nsnr) + (h - Hest(:, BICo)).' ...
          * (h - Hest(:, BICo));
      MSE_multiBIC(nsnr) = MSE_multiBIC(nsnr) + (h - hmultiBIC).' ...
          * (h - hmultiBIC);
      MSE_AIC(nsnr) = MSE_AIC(nsnr) + (h - Hest(:, AICo)).' ...
          * (h - Hest(:, AICo));
      MSE_AICc(nsnr) = MSE_AICc(nsnr) + (h - Hest(:, AICco)).' ...
          * (h - Hest(:, AICco));

      if L == (BICo)
         corrBIC(nsnr) = corrBIC(nsnr) + 1;
      end
      if L == (AICo)
         corrAIC(nsnr) = corrAIC(nsnr) + 1;
      end
      if L == (AICco)
         corrAICc(nsnr) = corrAICc(nsnr) + 1;
      end
      
      % BOSS/BPM for full a priori knowledge
      [h_mmse,h_oneterm,Lest] = bossbpm(X,y,s2(nsnr),diag(G), ...
                      ones(n,1)/n);
      
      MSE_Bayes(nsnr) = MSE_Bayes(nsnr) + (h - h_mmse).' * (h - ...
                                                        h_mmse);
      MSE_Bayes_oneterm(nsnr) = MSE_Bayes_oneterm(nsnr) + ...
          (h - h_oneterm).' * (h - h_oneterm);

      if L == Lest
        corrBayes(nsnr) = corrBayes(nsnr) + 1;
      end
      ordBayes(nsnr, Lest) = ordBayes(nsnr, Lest) + 1;
      if ((h-h_mmse)'*(h-h_mmse)) < min([(h-Hest(:,BICo))'*(h-Hest(:,BICo)), ...
                            (h-Hest(:,AICo))'*(h-Hest(:,AICo)), ...
                            (h-Hest(:,AICco))'*(h-Hest(:,AICco))])
         Bayes_beatsIC(nsnr) = ...
             Bayes_beatsIC(nsnr) + 1;
      end
   
      % HERE STARTS THE EMPIRICAL BAYES ESTIMATOR using Smoothed
      % estimator with true Psi matrix
        
      [hemp, hemp_oneterm,Lestemp] = bossbpm(X,y,s2hat_ub, ...
        rho, ones(n,1)/n);
      MSE_empBayesSmooth(nsnr) = MSE_empBayesSmooth(nsnr) + ...
          (h - hemp).' * (h - hemp);
      MSE_empBayesSmooth_oneterm(nsnr) = MSE_empBayesSmooth_oneterm(nsnr) + ...
          (h - hemp_oneterm).' * (h - hemp_oneterm);
      if L == Lestemp
        corrempBayesSmooth(nsnr) = corrempBayesSmooth(nsnr) + 1;
      end
      ordempBayesSmooth(nsnr, Lestemp) = ordempBayesSmooth(nsnr, Lestemp) + 1;
      if ((h-hemp)'*(h-hemp)) < min([(h-Hest(:,BICo))'*(h-Hest(:,BICo)), ...
                            (h-Hest(:,AICo))'*(h-Hest(:,AICo)), ...
                            (h-Hest(:,AICco))'*(h-Hest(:,AICco))])
         empBayesSmooth_beatsIC(nsnr) = ...
             empBayesSmooth_beatsIC(nsnr) + 1;
      end
      
      % HERE STARTS THE EMPIRICAL BAYES ESTIMATOR using Smoothed
      % estimator with Fourier Psi matrix
        
      [hemp, hemp_oneterm,Lestemp] = bossbpm(X,y,s2hat_ub, ...
        rhoF, ones(n,1)/n);
      MSE_empBayesSmoothF(nsnr) = MSE_empBayesSmoothF(nsnr) + ...
          (h - hemp).' * (h - hemp);
      MSE_empBayesSmoothF_oneterm(nsnr) = MSE_empBayesSmoothF_oneterm(nsnr) + ...
          (h - hemp_oneterm).' * (h - hemp_oneterm);
      if L == Lestemp
        corrempBayesSmoothF(nsnr) = corrempBayesSmoothF(nsnr) + 1;
      end
      ordempBayesSmoothF(nsnr, Lestemp) = ordempBayesSmoothF(nsnr, Lestemp) + 1;
      if ((h-hemp)'*(h-hemp)) < min([(h-Hest(:,BICo))'*(h-Hest(:,BICo)), ...
                            (h-Hest(:,AICo))'*(h-Hest(:,AICo)), ...
                            (h-Hest(:,AICco))'*(h-Hest(:,AICco))])
         empBayesSmoothF_beatsIC(nsnr) = ...
             empBayesSmoothF_beatsIC(nsnr) + 1;
      end


      % HERE STARTS THE EMPIRICAL BAYES ESTIMATOR using Smoothed
      % estimator with penalty on the Difference of rho
        
      [hemp, hemp_oneterm,Lestemp] = bossbpm(X,y,s2hat_ub, ...
        rhoD, ones(n,1)/n);
      MSE_empBayesSmoothD(nsnr) = MSE_empBayesSmoothD(nsnr) + ...
          (h - hemp).' * (h - hemp);
      MSE_empBayesSmoothD_oneterm(nsnr) = MSE_empBayesSmoothD_oneterm(nsnr) + ...
          (h - hemp_oneterm).' * (h - hemp_oneterm);
      if L == Lestemp
        corrempBayesSmoothD(nsnr) = corrempBayesSmoothD(nsnr) + 1;
      end
      ordempBayesSmoothD(nsnr, Lestemp) = ordempBayesSmoothD(nsnr, Lestemp) + 1;
      if ((h-hemp)'*(h-hemp)) < min([(h-Hest(:,BICo))'*(h-Hest(:,BICo)), ...
                            (h-Hest(:,AICo))'*(h-Hest(:,AICo)), ...
                            (h-Hest(:,AICco))'*(h-Hest(:,AICco))])
         empBayesSmoothD_beatsIC(nsnr) = ...
             empBayesSmoothD_beatsIC(nsnr) + 1;
      end
   
      % HERE STARTS THE EMPIRICAL BAYES ESTIMATOR using Smoothed
      % estimator with penalty on the Difference of rho obtained
      % using GCV
        
      [hemp, hemp_oneterm,Lestemp] = bossbpm(X,y,s2hat_ub, ...
        rhoGCV, ones(n,1)/n);
      MSE_empBayesSmoothGCV(nsnr) = MSE_empBayesSmoothGCV(nsnr) + ...
          (h - hemp).' * (h - hemp);
      MSE_empBayesSmoothGCV_oneterm(nsnr) = MSE_empBayesSmoothGCV_oneterm(nsnr) ...
          + (h - hemp_oneterm).' * (h - hemp_oneterm);
      if L == Lestemp
        corrempBayesSmoothGCV(nsnr) = corrempBayesSmoothGCV(nsnr) + 1;
      end
      ordempBayesSmoothGCV(nsnr, Lestemp) = ...
          ordempBayesSmoothGCV(nsnr, Lestemp) + 1;
      if ((h-hemp)'*(h-hemp)) < min([(h-Hest(:,BICo))'*(h-Hest(:,BICo)), ...
                            (h-Hest(:,AICo))'*(h-Hest(:,AICo)), ...
                            (h-Hest(:,AICco))'*(h-Hest(:,AICco))])
         empBayesSmoothGCV_beatsIC(nsnr) = ...
             empBayesSmoothGCV_beatsIC(nsnr) + 1;
      end
   
   end
end

total_time = toc

save(['Results_Nt', num2str(Nt), 'n', num2str(n)])

fig1 = figure(1);
clf
semilogy(s2_db, [MSE_AICc / Nmonte, MSE_BIC / Nmonte, ...
    MSE_multiBIC / Nmonte, MSE_Bayes / Nmonte, ...
                 MSE_empBayesSmooth / Nmonte, ...
                 MSE_empBayesSmoothF / Nmonte, ...
                 MSE_empBayesSmoothD / Nmonte])
hold on
semilogy(s2_db, MSE_empBayesSmoothGCV / Nmonte, 'k--o')
legend('LS with AICc-order', 'LS with BIC-order', 'BIC model avg', 'BPM', ... 
       'BPM emp, true \Psi', ...
       ['BPM emp, Fourier \Psi, r = ', num2str(sum(kF)+1)], ...
       ['BPM emp, \lambda = ', num2str(lam)], ...
       'BPM emp, GCV')
xlabel('noise variance \sigma^2 [dB]')
ylabel('MSE')
axis('tight')
ax1 = get(fig1,'CurrentAxes');
set(ax1,'Xdir','reverse')

fig2 = figure(2);
clf
plot(s2_db, [corrBIC / Nmonte, corrAICc / Nmonte, corrBayes / Nmonte, ...
    corrempBayesSmooth / Nmonte, corrempBayesSmoothF/Nmonte, ...
             corrempBayesSmoothD/Nmonte]*100)
hold on
plot(s2_db, [corrempBayesSmoothGCV/Nmonte]*100, 'k--o')
legend('BIC','AICc','BOSS', 'BOSS emp, true \Psi', ...
       ['BOSS emp Fourier \Psi, r = ', num2str(sum(kF)+1)], ...
       ['BOSS emp, \lambda = ', num2str(lam)], ...
       'BOSS emp, GCV')
xlabel('noise variance \sigma^2')
ylabel('percentage of correctly selected order')
ax2 = get(fig2,'CurrentAxes');
set(ax2,'Xdir','reverse')