www.pudn.com > trackingdemos.zip > MDL_demo2.m

% MDL_demo2.m 
% 
% demo version for comparison between the MDL model selection and the GLRT approach 
% in a one-dimensional detection problem using ML-PDA (cont. called by MDL_demo.m) 
disp('To compare MDL vs. GLRT, Monte Carlo simulations are needed'); 
disp('We consider three situations: '); 
disp('             1. No signal presents'); 
disp('             2. One signal presents at location 3'); 
disp('             3. Two signals present at locations -3 and 4'); 
disp('100 MC runs for each scenario and the numbers of making correct decisions'); 
disp('are listed ...'); 
disp('»'); 
disp('»'); 
disp('pause'); 
pause(5); 
disp('»'); 
disp('»'); 
disp('Number of scans = 40'); 
disp('Pd = 0.7'); 
disp('std of measurement noise = 1'); 
disp('Expected number of FAs per unit volume = 0.08'); 
disp('»'); 
disp('»'); 
disp('The simulation may take quite a while, please wait...'); 
lambda = 0.08; 
runs = 100; 
theta1 = 4; 
theta2 = -3; 
 
index1 = 0; 
index2 = 0; 
index3 = 0; 
index4 = 0; 
test_statH0 = []; 
test_statH1 = []; 
test_statH2 = []; 
test_newH0 = []; 
test_newH1 = []; 
test_newH2 = []; 
test_newnewH2 = []; 
h = waitbar(0, 'Please wait...'); 
for i=1:runs 
   for j=1:N 
      Z_H0(j).Z = []; 
      Z_H1(j).Z = []; 
      Z_H2(j).Z = []; 
      M1 = poissrnd(2*A*lambda); 
      M2 = poissrnd(2*A*lambda); 
      M3 = poissrnd(2*A*lambda); 
      if M1 > 0 
         Z_H0(j).Z = [Z_H0(j).Z, 2*A.*(rand(1,M1)-.5)]; 
      end 
      if M2 > 0 
         Z_H1(j).Z = [Z_H1(j).Z, 2*A.*(rand(1,M2)-.5)]; 
      end 
      if M3 > 0 
         Z_H2(j).Z = [Z_H2(j).Z, 2*A.*(rand(1,M3)-.5)]; 
      end 
       
      if rand < Pd 
         Z_H1(j).Z = [Z_H1(j).Z, theta + sigma.*randn]; 
      end 
      if rand < Pd 
         Z_H2(j).Z = [Z_H2(j).Z, theta1 + sigma.*randn]; 
      end 
      if rand < Pd  
         Z_H2(j).Z = [Z_H2(j).Z, theta2 + sigma.*randn]; 
      end 
   end 
         
   % test p(H2|H0), p(H1|H0), p(H0|H0) 
   theta_H0(i) = fminbnd('LLR_Gau_PDA2', -A, A, [], Z_H0, sigma, lambda, Pd, Pg); 
   NLLR_H0(i) = -LLR_Gau_PDA2(theta_H0(i), Z_H0, sigma, lambda, Pd, Pg); 
   number = get_meas_number(Z_H0); 
   test_statH0(i) = NLLR_H0(i) - .5*log(number/2/pi) - log(2*A/sigma); 
   if test_statH0(i) > 0 
      index1 = index1 + 1; 
      Z_H0_new = LLR_deduct_meas(Z_H0, theta_H0(i), sigma); 
      theta_newH0(index1) = fminbnd('LLR_Gau_PDA2', -A, A, [], Z_H0_new, sigma, lambda, Pd, Pg); 
      number = get_meas_number(Z_H0_new); 
      test_newH0(index1) = -LLR_Gau_PDA2(theta_newH0(index1), Z_H0_new, sigma, lambda, Pd, Pg) - .5*log(number/2/pi); 
   end 
        
   % test p(H2|H1), p(H1|H1), p(H0|H1) 
   theta_H1(i) = fminbnd('LLR_Gau_PDA2', -A, A, [], Z_H1, sigma, lambda, Pd, Pg); 
   NLLR_H1(i) = -LLR_Gau_PDA2(theta_H1(i), Z_H1, sigma, lambda, Pd, Pg); 
   number = get_meas_number(Z_H1); 
   test_statH1(i) = NLLR_H1(i) - .5*log(number/2/pi) - log(2*A/sigma); 
   if test_statH1(i) > 0 
      index2 = index2 + 1; 
      Z_H1_new = LLR_deduct_meas(Z_H1, theta_H1(i), sigma); 
      theta_newH1(index2) = fminbnd('LLR_Gau_PDA2', -A, A, [], Z_H1_new, sigma, lambda, Pd, Pg); 
      number = get_meas_number(Z_H1_new); 
      test_newH1(index2) = -LLR_Gau_PDA2(theta_newH1(index2), Z_H1_new, sigma, lambda, Pd, Pg) - .5*log(number/2/pi); 
   end 
 
	% test p(H3|H2), p(H2|H2), p(H1|H2), p(H0|H2) 
   theta_H2(i) = fminbnd('LLR_Gau_PDA2', -A, A, [], Z_H2, sigma, lambda, Pd, Pg); 
   NLLR_H2(i) = -LLR_Gau_PDA2(theta_H2(i), Z_H2, sigma, lambda, Pd, Pg); 
   number = get_meas_number(Z_H2); 
   test_statH2(i) = NLLR_H2(i) - .5*log(number/2/pi) - log(2*A/sigma); 
   if test_statH2(i) > 0 
      index3 = index3 + 1; 
      Z_H2_new = LLR_deduct_meas(Z_H2, theta_H2(i), sigma); 
      theta_newH2(index3) = fminbnd('LLR_Gau_PDA2', -A, A, [], Z_H2_new, sigma, lambda, Pd, Pg); 
      number = get_meas_number(Z_H2_new); 
      test_newH2(index3) = -LLR_Gau_PDA2(theta_newH2(index3), Z_H2_new, sigma, lambda, Pd, Pg) - .5*log(number/2/pi); 
      if test_newH2(index3) > 0 
         index4 = index4 + 1; 
         Z_H2_new_new = LLR_deduct_meas(Z_H2_new, theta_newH2(index3), sigma); 
         theta_newnewH2(index4) = fminbnd('LLR_Gau_PDA2', -A, A, [], Z_H2_new_new, sigma, lambda, Pd, Pg); 
         number = get_meas_number(Z_H2_new_new); 
         test_newnewH2(index4) = -LLR_Gau_PDA2(theta_newnewH2(index4), Z_H2_new_new, sigma, lambda, Pd, Pg) - .5*log(number/2/pi); 
      end 
   end 
   waitbar(i/runs); 
end 
close(h); 
P0_fct1 = length(find(test_statH0>0)); 
P0_fct2 = length(find(test_newH0>0)); 
P1_act1 = length(find(test_statH1>0)); 
P1_fct2 = length(find(test_newH1>0)); 
P2_fct1 = length(find(test_statH2>0)); 
P2_act2 = length(find(test_newH2>0)); 
P2_fct3 = length(find(test_newnewH2>0)); 
disp('Simulation results: '); 
disp(' '); 
disp('       Using MDL criterion'); 
disp('       1. No signal case'); 
string1 = ['          Number of correct decisions = ', num2str(runs-P0_fct1)]; 
disp(string1); 
disp('       2. One signal case'); 
string2 = ['          Number of correct decisions = ', num2str(P1_act1)]; 
disp(string2); 
disp('       3. Two signal case'); 
string3 = ['          Number of correct decisions = ', num2str(P2_act2)]; 
disp(string3); 
disp(' '); 
disp('       Using GLRT with 95% power under one signal case'); 
detection_threshold = prctile(NLLR_H1, 5); 
string4 = ['          Detection threshold for one signal vs. null: ', num2str(round(detection_threshold))]; 
disp(string4); 
disp(' '); 
disp('       1. No signal case'); 
string5 = ['          Number of correct decisions = ', num2str(length(find(NLLR_H0