From 5ac5ef1787a7100f1a7e724a780ebd10e4b6cf6c Mon Sep 17 00:00:00 2001 From: Ksyer <> Date: Tue, 16 Jul 2024 15:56:53 +0800 Subject: [PATCH] Update FAR vs PD --- FAR vs PD/FAR_simu.m | 122 ++++++++++++++++++++++++++++++----------- FAR vs PD/Main.m | 15 ++--- FAR vs PD/check_orth.m | 25 +++++++++ FAR vs PD/recovery.m | 10 +++- 4 files changed, 133 insertions(+), 39 deletions(-) create mode 100755 FAR vs PD/check_orth.m diff --git a/FAR vs PD/FAR_simu.m b/FAR vs PD/FAR_simu.m index c62a710..0b2fc5d 100755 --- a/FAR vs PD/FAR_simu.m +++ b/FAR vs PD/FAR_simu.m @@ -1,4 +1,4 @@ -function P_fa = FAR_simu(threshold) +function threshold = FAR_simu(P_fa_FAR) global N M epi trail_times % Define measurement matrix @@ -9,64 +9,124 @@ function P_fa = FAR_simu(threshold) Lambda_C = setdiff(1:N*M, Lambda); % Try "trail_times" times - NN = N; - results = zeros(trail_times, N * M); - - if exist("FAR_recovery_1_LASSO_5.mat") - load("FAR_recovery_0_debiasedLASSO_500.mat", "results"); + results_0 = zeros(trail_times, N * M); + results_1 = zeros(trail_times, N * M); + thresholds_0 = []; + thresholds_1 = []; + + noise_sigma = 0.01; + noise_sigma_2 = noise_sigma ^ 2; + + if exist("520BP.mat") + load("50.mat", "results_0", "results_1"); else + % 0 假设 for T = 1: trail_times % y = Ax + n - n = randn(NN, 1) * 0.1; + n = randn(N, 1) * noise_sigma; + y_noise = n; + [x_hat, threshold] = debiased_LASSO(A, y_noise, P_fa_FAR, noise_sigma_2); + results_0(T, :) = x_hat; + thresholds_0 = [thresholds_0 threshold]; + end + % 1 假设 + for T = 1: trail_times + % y = Ax + n + n = randn(N, 1) * noise_sigma; y_noise = A * x + n; - - % x_hat = CS(A, y) - x_hat = recovery(A, y_noise); - results(T, :) = x_hat; + [x_hat, threshold] = debiased_LASSO(A, y_noise, P_fa_FAR, noise_sigma_2); + results_1(T, :) = x_hat; + thresholds_1 = [thresholds_1 threshold]; end end % Distribute of H_0 and H_1 - H_0_distribute = zeros(1, length(Lambda_C)); - H_1_distribute = zeros(1, length(Lambda)); + H_00_distribute = zeros(1, length(Lambda_C)); + H_01_distribute = zeros(1, length(Lambda)); + H_10_distribute = zeros(1, length(Lambda_C)); + H_11_distribute = zeros(1, length(Lambda)); i1 = 1; i2 = 1; - Ts = []; + T_00 = []; + T_01 = []; + T_10 = []; + T_11 = []; - for t = 1: trail_times - x_hat = results(t, :); - T = 0; + for threshold = 1: trail_times + x_0_hat = results_0(threshold, :); + x_1_hat = results_1(threshold, :); + t_00 = 0; + t_01 = 0; + t_10 = 0; + t_11 = 0; + for i = 1: length(x) if ismember(i, Lambda_C) - H_0_distribute(i1) = x_hat(i); + t_00 = t_00 + abs(x_0_hat(i)); + t_10 = t_10 + abs(x_1_hat(i)); + H_00_distribute(i1) = x_0_hat(i); + H_10_distribute(i1) = x_1_hat(i); i1 = i1 + 1; elseif ismember(i, Lambda) - T = T + abs(x_hat(i)); - H_1_distribute(i2) = x_hat(i); + t_01 = t_01 + abs(x_0_hat(i)); + t_11 = t_11 + abs(x_1_hat(i)); + H_01_distribute(i2) = x_0_hat(i); + H_11_distribute(i2) = x_1_hat(i); i2 = i2 + 1; end end - Ts = [Ts T]; + T_00 = [T_00 t_00]; + T_01 = [T_01 t_01]; + T_10 = [T_10 t_10]; + T_11 = [T_11 t_11]; end % Draw figure(1) - subplot(2, 1, 1); + subplot(2, 2, 1); title("Freq histogram of H_0"); - histfit(real(H_0_distribute)); + histfit(real(H_00_distribute)); - subplot(2, 1, 2); + subplot(2, 2, 2); title("Freq histogram of H_1"); - histfit(real(H_1_distribute)); + histfit(real(H_01_distribute)); + + subplot(2, 2, 3); + title("Freq histogram of H_0"); + histfit(real(H_10_distribute)); - H_0_mean = mean(H_0_distribute); - H_0_std = std(H_0_distribute); - H_1_mean = mean(H_1_distribute); - H_1_std = std(H_1_distribute); + subplot(2, 2, 4); + title("Freq histogram of H_1"); + histfit(real(H_11_distribute)); + + H_00_mean = mean(H_00_distribute); + H_00_std = std(H_00_distribute); + H_01_mean = mean(H_01_distribute); + H_01_std = std(H_01_distribute); + H_10_mean = mean(H_10_distribute); + H_10_std = std(H_10_distribute); + H_11_mean = mean(H_11_distribute); + H_11_std = std(H_11_distribute); + + fprintf("H_00: mu = %f, std = %f\n", H_00_mean, H_00_std); + fprintf("H_01: mu = %f, std = %f\n", H_01_mean, H_01_std); + fprintf("H_10: mu = %f, std = %f\n", H_10_mean, H_10_std); + fprintf("H_11: mu = %f, std = %f\n", H_11_mean, H_11_std); + + figure(2); + subplot(2, 2, 1); + histfit(real(T_00)); + + subplot(2, 2, 2); + histfit(real(T_01)); + + subplot(2, 2, 3); + histfit(real(T_10)); + + subplot(2, 2, 4); + histfit(real(T_11)); - fprintf("H_0: mu = %f, std = %f\n", H_0_mean, H_0_std); - fprintf("H_1: mu = %f, std = %f\n", H_1_mean, H_1_std); % P_fa = normcdf(threshold, H_0_mean, H_0_std); end diff --git a/FAR vs PD/Main.m b/FAR vs PD/Main.m index c58632e..37ac3ed 100755 --- a/FAR vs PD/Main.m +++ b/FAR vs PD/Main.m @@ -7,19 +7,20 @@ N = 64; M = 4; K = 10; lambda = zeros(M * N, 1); -lambda(:) = 0.3; -tau = 1e-4; -iter_max = 500; -% method = "debiased_LASSO"; +lambda(:) = 0.8; +tau = 1e-6; +iter_max = 50; +method = "debiased_LASSO"; % method = "LASSO"; -method = "BP"; +% method = "BP"; % method = "VAMP"; -trail_times = 1e2; +trail_times = 1000; epi = 0; extend_target = 1; %% FAR -P_fa_FAR = FAR_simu(3); +P_fa_FAR = 1e-5; +threshold_FAR = FAR_simu(P_fa_FAR); % fprintf("FAR: %f\n", P_fa_FAR); %% PD diff --git a/FAR vs PD/check_orth.m b/FAR vs PD/check_orth.m new file mode 100755 index 0000000..07216f7 --- /dev/null +++ b/FAR vs PD/check_orth.m @@ -0,0 +1,25 @@ +clc; clear; +N = 16; +M = 4; +epi = 0; + +d_n = floor(rand(N, 1)*M) / M; +A = get_Psi(N, M, d_n, epi); +B = get_Psi_2(N, M, d_n, epi); +C = A - B; +real_C = real(C); +imag_C = imag(C); + +% subplot(121) +% heatmap(real_C) +% +% subplot(122) +% heatmap(imag_C) + +r = zeros(N); +for i = 1: N + for j = 1: N + r(i, j) = A(i, :) * (A(j, :)'); + end +end +heatmap(abs(r)); diff --git a/FAR vs PD/recovery.m b/FAR vs PD/recovery.m index 3132ca0..a69d26b 100755 --- a/FAR vs PD/recovery.m +++ b/FAR vs PD/recovery.m @@ -25,6 +25,14 @@ function x_hat = recovery(A, y_noise) x_d_CROD = x_LASSO + A'*(y_noise - A*x_LASSO)/Q_hat; x_hat = x_d_CROD; + RSS = 1/sz(2) * norm(y_noise - A * x_LASSO, 2)^2; + noise_sigma_2 = 1e-4; + P_fa = normcdf(1, 0, 1); + sigma_w_2 = (gamma * (1-gamma)) / ((gamma - Rho)^2) * RSS + noise_sigma_2; + k_d = -sigma_w_2 * log(P_fa); + + fprintf("%f %f %f\n", gamma, sigma_w_2, k_d); + elseif method == "LASSO" sz = size(A); N = sz(2); @@ -47,7 +55,7 @@ function x_hat = recovery(A, y_noise) cvx_end else global lambda tau iter_max; - [x_hat, z_hat_d] = cVAMPro(y_noise, A, lambda, tau, iter_max); + [x_hat, x_hat_d] = cVAMPro(y_noise, A, lambda, tau, iter_max); end end