From 9cbe8e31eacaebdcc7e8023427e2cb3bc0ccb6f6 Mon Sep 17 00:00:00 2001 From: Ksyer <> Date: Mon, 22 Jul 2024 21:17:57 +0800 Subject: [PATCH] Update wide vs narrow code --- wide vs narrow/Main3.m | 253 ++++++++++++++++ .../single_point_test/All_Recovery.m | 104 +++++++ .../single_point_test/All_Recovery2.m | 95 ++++++ wide vs narrow/single_point_test/get_lambda.m | 19 ++ .../single_point_test/get_sigma_var.m | 71 +++++ .../single_point_test/get_sparse_vector.m | 28 ++ wide vs narrow/single_point_test/parameters.m | 88 ++++++ wide vs narrow/single_point_test/query2.m | 39 +++ wide vs narrow/single_point_test/recovery.m | 52 ++++ .../signle_point_mento_carlo_test.m | 269 +++++++++++++++++ .../signle_point_mento_carlo_train.m | 53 ++++ wide vs narrow/single_point_test/simulator.m | 196 +++++++++++++ wide vs narrow/single_point_test/swtest.m | 273 ++++++++++++++++++ .../single_point_test/test_lambda.m | 102 +++++++ 14 files changed, 1642 insertions(+) create mode 100755 wide vs narrow/Main3.m create mode 100755 wide vs narrow/single_point_test/All_Recovery.m create mode 100755 wide vs narrow/single_point_test/All_Recovery2.m create mode 100755 wide vs narrow/single_point_test/get_lambda.m create mode 100755 wide vs narrow/single_point_test/get_sigma_var.m create mode 100755 wide vs narrow/single_point_test/get_sparse_vector.m create mode 100755 wide vs narrow/single_point_test/parameters.m create mode 100755 wide vs narrow/single_point_test/query2.m create mode 100755 wide vs narrow/single_point_test/recovery.m create mode 100755 wide vs narrow/single_point_test/signle_point_mento_carlo_test.m create mode 100755 wide vs narrow/single_point_test/signle_point_mento_carlo_train.m create mode 100644 wide vs narrow/single_point_test/simulator.m create mode 100755 wide vs narrow/single_point_test/swtest.m create mode 100755 wide vs narrow/single_point_test/test_lambda.m diff --git a/wide vs narrow/Main3.m b/wide vs narrow/Main3.m new file mode 100755 index 0000000..b94e4b6 --- /dev/null +++ b/wide vs narrow/Main3.m @@ -0,0 +1,253 @@ +clc; clear; close all; + +parameters; + +P_fas = [1e-7, 5e-7, 1e-6, 5e-6, 1e-5, 5e-5, 1e-4, 5e-4, 1e-3, 5e-3, 1e-2, 5e-2, 1e-1, 5e-1, 1e0]; + +P_d_FARs = zeros(length(P_fas), 1); +P_d_narrows = zeros(length(P_fas), 1); + + +%% FAR + +% Signal Model +% epi = B / f_c; +epi = 0; +[C_n, A] = get_Psi(FAR_N, FAR_M, epi); +fac = A * A'; +A = A / sqrt(abs(fac(1, 1))); +f_n = f_c + C_n * B / M; + +[Lambda, Lambda_C, x] = get_sparse_vector(FAR_N, FAR_M, betas_wide, true); + +% Define filename +filename = ... + "data2/" + ... + "FAR" + "_" + ... + string(FAR_N) + "_" + ... + string(FAR_M)+ "_" + ... + method + "_" + ... + string(noise_sigma) + "_" + ... + string(length(Lambda)) + "_" + ... + string(LASSO_lambda) + "_" + ... + trail_times + ... + ".mat"; + + +for P_fa_idx = 1: length(P_fas) + P_fa = P_fas(P_fa_idx); + if exist(filename, "file") + load( ... + filename, ... + "recovery_results", "sigma_w2", "thresholds",... + "C_n", "A", "Lambda", "Lambda_C", "x" ... + ); + else + % Mento Carlo Recovery + [recovery_results, sigma_w2, thresholds] = All_Recovery(A, x, P_fa, noise_sigma / FAR_M); + + save( ... + filename, ... + "recovery_results", "sigma_w2", "thresholds",... + "C_n", "A", "Lambda", "Lambda_C", "x" ... + ); + end + + % Compare sigmas (from paper) and vars (from simulation) + % [ss, vs, H0_mean, H0_max, H1_mean, H1_max] = get_sigma_var(recovery_results, sigma_w2, x); + + T_opt = zeros(2, 1); + + delta_r = c / (2 * B); + tmp_1 = repelem(f_n, FAR_M); + tmp_2 = repmat(1: FAR_M, 1, FAR_N)'; + tmp_3 = exp(-1 * 1j * 2 * pi * 2 * delta_r .* tmp_2 .* tmp_1 / c); + + idx_1 = 1; + idx_2 = 1; + i1 = 1; i2 = 1; + + Lambda_distributes = zeros(2, length(Lambda) * trail_times); + Lambda_C_distributes = zeros(2, length(Lambda_C) * trail_times); + + for t = 1:trail_times + % 0 假设 + x_0_hat = squeeze(recovery_results(1, t, :)); + s = x_0_hat .* tmp_3; + for i = 1: FAR_N + L = (i - 1) * FAR_M + 1; + R = L - 1 + FAR_M; + T_opt(1, idx_1) = foo(s(L:R)); + idx_1 = idx_1 + 1; + end + + % 1 假设 + x_1_hat = squeeze(recovery_results(2, t, :)); + s = x_1_hat .* tmp_3; + T_opt(2, idx_2) = foo(s(Lambda)); + idx_2 = idx_2 + 1; + + for i = 1: length(Lambda) + length(Lambda_C) + if ismember(i, Lambda_C) + Lambda_C_distributes(1, i1) = x_0_hat(i); + Lambda_C_distributes(2, i1) = x_1_hat(i); + i1 = i1 + 1; + elseif ismember(i, Lambda) + Lambda_distributes(1, i2) = x_0_hat(i); + Lambda_distributes(2, i2) = x_1_hat(i); + i2 = i2 + 1; + end + end + end + + % figure; + % subplot(211); plot(real(sum(squeeze(recovery_results(1, 1:10, :))))); + % subplot(212); plot(real(sum(squeeze(recovery_results(2, 1:10, :))))); + + % figure; + % subplot(211); histfit(T_opt(1, 1:idx_1 - 1)); + % subplot(212); histfit(T_opt(2, 1:idx_2 - 1)); + + figure; + subplot(2, 2, 1); histfit(real(Lambda_C_distributes(1, :))); + subplot(2, 2, 2); histfit(real(Lambda_distributes(1, :))); + subplot(2, 2, 3); histfit(real(Lambda_C_distributes(2, :))); + subplot(2, 2, 4); histfit(real(Lambda_distributes(2, :))); + + figure; + histfit(real([Lambda_distributes(2, :) Lambda_C_distributes(2, :)])); + + means = [ ... + mean(Lambda_C_distributes(1, :)), ... + mean(Lambda_distributes(1, :)), ... + mean(Lambda_C_distributes(2, :)), ... + mean(Lambda_distributes(2, :)), ... + mean([Lambda_distributes(1, :) Lambda_C_distributes(1, :)]), ... + mean([Lambda_distributes(1, :) Lambda_C_distributes(1, :) Lambda_C_distributes(2, :)]), ... + mean([Lambda_distributes(2, :)]) ... + ]; + + stds = [ ... + std(Lambda_C_distributes(1, :)), ... + std(Lambda_distributes(1, :)), ... + std(Lambda_C_distributes(2, :)), ... + std(Lambda_distributes(2, :)), ... + std([Lambda_distributes(1, :) Lambda_C_distributes(1, :)]), ... + std([Lambda_distributes(1, :) Lambda_C_distributes(1, :) Lambda_C_distributes(2, :)]), ... + std([Lambda_distributes(2, :)]) ... + ]; + + + fprintf("H_00: mu = %.15f, std = %.15f\n", means(1), stds(1)); + fprintf("H_01: mu = %.15f, std = %.15f\n", means(2), stds(2)); + fprintf("H_10: mu = %.15f, std = %.15f\n", means(3), stds(3)); + fprintf("H_11: mu = %.15f, std = %.15f\n", means(4), stds(4)); + fprintf("H_0: mu = %.15f, std = %.15f\n", means(5), stds(5)); + fprintf("H_0 + Lambda^C: mu = %.15f, std = %.15f\n", means(6), stds(6)); + fprintf("H_1 Lambda: mu = %.15f, std = %.15f\n\n", means(7), stds(7)); + + + H_0_mean = mean(T_opt(1, 1:idx_1-1)); + H_0_std = std(T_opt(1, 1:idx_1-1)); + H_1_mean = mean(T_opt(2, 1:idx_2-1)); + H_1_std = std(T_opt(2, 1:idx_2-1)); + + + threshold = norminv(1 - P_fa, H_0_mean, H_0_std); + P_d = 1 - normcdf(threshold, H_1_mean, H_1_std); + P_d_FARs(P_fa_idx) = P_d; +end + +format long; +fprintf("H_0 mean: %.15f\n", H_0_mean); +fprintf("H_0 std: %.15f\n", H_0_std); +fprintf("H_1 mean: %.15f\n", H_1_mean); +fprintf("H_1 std: %.15f\n", H_1_std); + +%% PD +[s_T_narrow, A] = narrow_signal_model_2(FAR_N); +[Lambda, Lambda_C, x] = get_sparse_vector(FAR_N, 1, betas_narrow, false); +filename = ... + "data2/" + ... + "PD" + "_" + ... + method + "_" + ... + string(noise_sigma) + "_" + ... + string(length(Lambda)) + "_" + ... + trail_times + ... + ".mat"; + +for P_fa_idx = 1: length(P_fas) + [s_T, A] = narrow_signal_model_2(FAR_N); + P_fa = P_fas(P_fa_idx); + if exist(filename, "file") + load( ... + filename, ... + "recovery_results", "sigma_w2", "thresholds",... + "f_c", "A", "Lambda", "Lambda_C", "x" ... + ); + else + % Mento Carlo Recovery + [recovery_results, sigma_w2, thresholds] = All_Recovery(A, x, P_fa, noise_sigma); + + save( ... + filename, ... + "recovery_results", "sigma_w2", "thresholds",... + "f_c", "A", "Lambda", "Lambda_C", "x" ... + ); + end + + T_opt = zeros(2, trail_times); + + delta_r = T_p * c / 2; + tmp_1 = (1: FAR_N)'; + tmp_3 = exp(-1 * 1j * 2 * pi * 2 * delta_r * f_c .* tmp_1 / c); + + for t = 1:trail_times + % 0 假设 + x_hat = recovery_results(1, t, :); + s = squeeze(x_hat) .* tmp_3; + for i = 1: FAR_N + T_opt(1, idx_1) = foo(s(i)); + idx_1 = idx_1 + 1; + end + + % 1 假设 + x_hat = recovery_results(2, t, :); + s = squeeze(x_hat) .* tmp_3; + T_opt(2, idx_2) = foo(s(Lambda)); + idx_2 = idx_2 + 1; + end + + % figure; + % subplot(211); plot(real(sum(squeeze(recovery_results(1, 1:10, :))))); + % subplot(212); plot(real(sum(squeeze(recovery_results(2, 1:10, :))))); + + H_0_mean = mean(T_opt(1, 1:idx_1-1)); + H_0_std = std(T_opt(1, 1:idx_1-1)); + H_1_mean = mean(T_opt(2, 1:idx_2-1)); + H_1_std = std(T_opt(2, 1:idx_2-1)); + + threshold = norminv(1 - P_fa, H_0_mean, H_0_std); + P_d = 1 - normcdf(threshold, H_1_mean, H_1_std); + P_d_narrows(P_fa_idx) = P_d; +end + +format long; +fprintf("H_0 mean: %.15f\n", H_0_mean); +fprintf("H_0 std: %.15f\n", H_0_std); +fprintf("H_1 mean: %.15f\n", H_1_mean); +fprintf("H_1 std: %.15f\n", H_1_std); + + +%% ROC + +figure; +hold on; +grid on; +semilogx(P_fas, P_d_FARs); +semilogx(P_fas, P_d_narrows); +xlim([8e-8, 1]); +% ylim([0, 1]); +xlabel("P_{fa}"); +ylabel("P_{d}"); +title("ROC"); diff --git a/wide vs narrow/single_point_test/All_Recovery.m b/wide vs narrow/single_point_test/All_Recovery.m new file mode 100755 index 0000000..1ffb26e --- /dev/null +++ b/wide vs narrow/single_point_test/All_Recovery.m @@ -0,0 +1,104 @@ +function [recovery_results, sigma_w2, thresholds] = All_Recovery(A, x, P_fa, noise_sigma) + + global trail_times method LASSO_lambda tau iter_max; + sz = size(A); + M = sz(1); + N = sz(2); + + recovery_results = zeros(2, trail_times, N); + sigma_w2 = zeros(2, trail_times, 1); + thresholds = zeros(2, trail_times, 1); + figure; + + for hypo = 2: -1: 1 + % hypo- + hypo = 3 - hypo; + h = waitbar(0, 'ڷ' + string(hypo-1) + ''); + for T = 1:trail_times + waitbar(T / trail_times, h); + + noise = get_noise(noise_sigma, M, 1); + + % y = Ax + n + if hypo == 1 + y_noise = noise; + else + y_noise = A * x + noise; + end + + if method == "debiased_LASSO" + [x_hat, sigma_w_2, threshold] = debiased_LASSO(A, y_noise, P_fa, noise_sigma^2, LASSO_lambda); + sigma_w2(hypo, T) = sigma_w_2; + thresholds(hypo, T) = threshold; + elseif method == "debiased_LASSO_FISTA" + [x_hat, sigma_w_2, threshold] = debiased_LASSO_FISTA(A, y_noise, P_fa, noise_sigma^2, LASSO_lambda); + sigma_w2(hypo, T) = sigma_w_2; + thresholds(hypo, T) = threshold; + elseif method == "cVAMPro" + [x_LASSO, x_hat_d] = cVAMPro(y_noise, A, LASSO_lambda, tau, 100); + % x_LASSO = FISTA(y_noise, A, LASSO_lambda, 1e-5); + + + sz = size(A); + n = sz(2); + gamma = sz(1) / sz(2); + lambda = LASSO_lambda; + + % CRODȥƫ + rho_active = sum(abs(x_LASSO) > 1e-3)/n; + Q_hat = (gamma - rho_active)/(1 - rho_active); + Rho = sum((abs(x_LASSO) > 1e-3).* (2 - lambda./(Q_hat*abs(x_LASSO) + lambda))) / 2 / n; + diff = 1; + while(diff > 1e-4) + Rho_pre = Rho; + Rho = sum((abs(x_LASSO) > 1e-3).* (2 - lambda./((gamma-Rho)/(1-Rho)*abs(x_LASSO) + lambda))) / 2 / n; + diff = abs(Rho - Rho_pre); + end + Q_hat = (gamma-Rho)/(1-Rho); + x_d_CROD = x_LASSO + A'*(y_noise - A*x_LASSO)/Q_hat; % x_d_CROD == x_hat_d + sigma_n = noise_sigma; + + % CROD޺ͼͳ + RSS = sum(abs(y_noise - A * x_LASSO).^2)/length(y_noise); + chi = Rho*(1 - Rho)/(gamma - Rho); + if chi ~= 0 + chi_temp = sqrt((chi+1)*(chi+1)-4*gamma*chi); + z = -(1 - chi + chi_temp) / (2*chi); + z_prime = -(1 - 2*gamma*chi + chi + chi_temp) / (2*chi*chi*chi_temp); + G_prime = (z + 1/chi); + G_wprime = (z_prime + 1/chi/chi); + chi_hat = gamma/2*G_wprime*RSS/(G_prime - chi*G_wprime)... + + (G_prime*G_prime/2 - gamma/2*G_wprime)*sigma_n*sigma_n/(G_prime - chi*G_wprime); + else + G_prime = gamma; + G_wprime = gamma*(1-gamma); + chi_hat = gamma/2*G_wprime*RSS/(G_prime - chi*G_wprime)... + + (G_prime*G_prime/2 - gamma/2*G_wprime)*sigma_n*sigma_n/(G_prime - chi*G_wprime); + end + sigma_CROD = sqrt(2*chi_hat) / Q_hat; + + sigma_w2(hypo, T) = sigma_CROD^2; + thresholds(hypo, T) = -sigma_CROD^2 * log(P_fa); + % figure; subplot(221); plot(real(x)); subplot(222); plot(real(x_LASSO)); subplot(223); plot(real(x_hat_d)); subplot(224); plot(real(x_d_CROD)); + if hypo == 2 + x_hat = x_hat_d; + % pp = real(x_hat - x); + % [is_not_norm, tmp, tmp] = swtest(pp, 0.1); + % + % if is_not_norm + % subplot(211); plot(pp); subplot(212); histfit(pp); + % is_not_norm + % end + else + x_hat = x_d_CROD; + end + else + x_hat = recovery(A, y_noise, method); + end + + + recovery_results(hypo, T, :) = x_hat; + end + close(h); + end +end diff --git a/wide vs narrow/single_point_test/All_Recovery2.m b/wide vs narrow/single_point_test/All_Recovery2.m new file mode 100755 index 0000000..7f5444b --- /dev/null +++ b/wide vs narrow/single_point_test/All_Recovery2.m @@ -0,0 +1,95 @@ +function [recovery_results, sigma_w2, thresholds] = All_Recovery2(A, x, P_fa, noise_sigma, trail_times, LASSO_lambda, tau) + + sz = size(A); + M = sz(1); + N = sz(2); + + recovery_results = zeros(2, trail_times, N); + sigma_w2 = zeros(2, trail_times, 1); + thresholds = zeros(2, trail_times, 1); + + for hypo = 2: -1: 1 + % hypo- + hypo = 3 - hypo; + h = waitbar(0, 'ڷ' + string(hypo-1) + ''); + for T = 1:trail_times + waitbar(T / trail_times, h); + + noise = get_noise(noise_sigma, M, 1); + + % y = Ax + n + if hypo == 1 + y_noise = noise; + else + y_noise = A * x + noise; + end + + [x_LASSO, x_hat_d] = cVAMPro(y_noise, A, LASSO_lambda, tau, 100); + % x_LASSO = FISTA(y_noise, A, LASSO_lambda, 1e-5); + + + sz = size(A); + n = sz(2); + gamma = sz(1) / sz(2); + lambda = LASSO_lambda; + + % CRODȥƫ + rho_active = sum(abs(x_LASSO) > 1e-3)/n; + Q_hat = (gamma - rho_active)/(1 - rho_active); + Rho = sum((abs(x_LASSO) > 1e-3).* (2 - lambda./(Q_hat*abs(x_LASSO) + lambda))) / 2 / n; + diff = 1; + while(diff > 1e-4) + Rho_pre = Rho; + Rho = sum((abs(x_LASSO) > 1e-3).* (2 - lambda./((gamma-Rho)/(1-Rho)*abs(x_LASSO) + lambda))) / 2 / n; + diff = abs(Rho - Rho_pre); + end + Q_hat = (gamma-Rho)/(1-Rho); + x_d_CROD = x_LASSO + A'*(y_noise - A*x_LASSO)/Q_hat; % x_d_CROD == x_hat_d + sigma_n = noise_sigma; + + % CROD޺ͼͳ + RSS = sum(abs(y_noise - A * x_LASSO).^2)/length(y_noise); + chi = Rho*(1 - Rho)/(gamma - Rho); + if chi ~= 0 + chi_temp = sqrt((chi+1)*(chi+1)-4*gamma*chi); + z = -(1 - chi + chi_temp) / (2*chi); + z_prime = -(1 - 2*gamma*chi + chi + chi_temp) / (2*chi*chi*chi_temp); + G_prime = (z + 1/chi); + G_wprime = (z_prime + 1/chi/chi); + chi_hat = gamma/2*G_wprime*RSS/(G_prime - chi*G_wprime)... + + (G_prime*G_prime/2 - gamma/2*G_wprime)*sigma_n*sigma_n/(G_prime - chi*G_wprime); + else + G_prime = gamma; + G_wprime = gamma*(1-gamma); + chi_hat = gamma/2*G_wprime*RSS/(G_prime - chi*G_wprime)... + + (G_prime*G_prime/2 - gamma/2*G_wprime)*sigma_n*sigma_n/(G_prime - chi*G_wprime); + end + sigma_CROD = sqrt(2*chi_hat) / Q_hat; + + if sigma_CROD > 1 + sigma_CROD = sigma_CROD; + end + sigma_w2(hypo, T) = sigma_CROD^2; + thresholds(hypo, T) = -sigma_CROD^2 * log(P_fa); + % figure; subplot(221); plot(real(x)); subplot(222); plot(real(x_LASSO)); subplot(223); plot(real(x_hat_d)); subplot(224); plot(real(x_d_CROD)); + + % % data4 folder + % if hypo == 2 + % x_hat = x_d_CROD; + % % x_hat = x_hat_d; + % else + % x_hat = x_d_CROD; + % end + + % data5 folder + if hypo == 2 + x_hat = x_hat_d; + else + x_hat = x_d_CROD; + end + + recovery_results(hypo, T, :) = x_hat; + end + close(h); + end +end diff --git a/wide vs narrow/single_point_test/get_lambda.m b/wide vs narrow/single_point_test/get_lambda.m new file mode 100755 index 0000000..eda8639 --- /dev/null +++ b/wide vs narrow/single_point_test/get_lambda.m @@ -0,0 +1,19 @@ +function [LASSO_lambda, recovery_results] = get_lambda(A, x, FAR_N, FAR_M, sigma, lambdas) + + noise = get_noise(sigma, FAR_N, 1); + y = A * x; + y_noise = y + noise; + + tau = 1e-6; + iter_max = 1000; + + recovery_results = zeros(length(lambdas), 1); + for lambda_idx = 1: length(lambdas) + LASSO_lambda = lambdas(lambda_idx); + [x_LASSO, x_hat_d] = cVAMPro(y_noise, A, LASSO_lambda, tau, iter_max); + recovery_results(lambda_idx) = sum(real(x_LASSO(FAR_M + 1: 2 * FAR_M))); + % if abs(recovery_results(lambda_idx) - FAR_M) / FAR_M < 0.01 + % return + % end + end +end \ No newline at end of file diff --git a/wide vs narrow/single_point_test/get_sigma_var.m b/wide vs narrow/single_point_test/get_sigma_var.m new file mode 100755 index 0000000..91323df --- /dev/null +++ b/wide vs narrow/single_point_test/get_sigma_var.m @@ -0,0 +1,71 @@ +function [sigmas, vars, H0_mean, H0_max, H1_mean, H1_max] = get_sigma_var(recovery_results, sigma_w2, x) + + sz = size(recovery_results); + trail_times = sz(2); + + sigmas = zeros(2, trail_times); + vars = zeros(2, trail_times); + id = 1; + + for T = 1:trail_times + + x_0_hat = squeeze(recovery_results(1, T, :)); + x_1_hat = squeeze(recovery_results(2, T, :)) - x; + + if anynan(x_1_hat) + continue; + end + + % if sigma_w2(2, T) > 0.1 + % continue; + % end + + sigmas(1, id) = sigma_w2(1, T); + sigmas(2, id) = sigma_w2(2, T); + + vars(1, id) = var(x_0_hat); + vars(2, id) = var(x_1_hat); + + id = id + 1; + end + sigmas(:, end - trail_times + id: end) = []; + vars(:, end - trail_times + id: end) = []; + + REE = abs(vars - sigmas) ./ vars; + H0_mean = mean(REE(1, :)); + H0_max = max(REE(1, :)); + H1_mean = mean(REE(2, :)); + H1_max = max(REE(2, :)); + + if 1 == 9 + figure; + subplot(211); + plot(sigmas(1, :)); + hold on; + plot(vars(1, :)); + xlabel("Mento Carlo Times"); + ylabel("\sigma^2") + legend("\sigma^2 (CROD)", "\sigma_w^2 (Expr)"); + subplot(212); + plot(sigmas(2, :)); + hold on; + plot(vars(2, :)); + xlabel("Mento Carlo Times"); + ylabel("\sigma^2"); + legend("\sigma^2 (CROD)", "\sigma_w^2 (Expr)"); + + figure; + subplot(211); + plot(REE(1, :)); + yline(H0_mean); + xlabel("Mento Carlo Times"); + ylabel("REE"); + legend("mean = " + string(H0_mean), "max = " + string(H0_max)); + subplot(212); + plot(REE(2, :)); + yline(H1_mean); + xlabel("Mento Carlo Times"); + ylabel("REE"); + legend("mean = " + string(H1_mean), "max = " + string(H1_max)); + end +end diff --git a/wide vs narrow/single_point_test/get_sparse_vector.m b/wide vs narrow/single_point_test/get_sparse_vector.m new file mode 100755 index 0000000..30d16d1 --- /dev/null +++ b/wide vs narrow/single_point_test/get_sparse_vector.m @@ -0,0 +1,28 @@ +function [Lambda, Lambda_C, x] = get_sparse_vector(N, M, Amp, extend_target) + x = zeros(N * M, 1); + + if nargin < 4 + extend_target = true; + end + + % Support Set + if length(x) > 2 * M && extend_target + Lambda = M + 1:2 * M; + elseif length(x) < 2 * M && extend_target + Lambda = 1:M; + else + Lambda = 2; + end + + if length(Amp) == 1 + x(Lambda) = Amp; + elseif length(Amp) <= max(Lambda) + x(Lambda) = Amp(Lambda); + else + % Warning + x(Lambda) = Amp(end); + end + + Lambda_C = setdiff(1:N*M, Lambda); + +end diff --git a/wide vs narrow/single_point_test/parameters.m b/wide vs narrow/single_point_test/parameters.m new file mode 100755 index 0000000..3346e40 --- /dev/null +++ b/wide vs narrow/single_point_test/parameters.m @@ -0,0 +1,88 @@ +global ... + B ... + T_p ... + f_s ... + T_s ... + f_c ... + PRF ... + T_r ... + t_single_pulse ... + c ... + N ... + M ... + R_0... + delta_R_wide ... + delta_R_narrow ... + N_wide ... + N_narrow ... + ranges_wide... + ranges_narrow ... + betas_wide ... + betas_narrow ... + trail_times... + method ... + N_high ... + LASSO_lambda ... + tau ... + iter_max ... + FAR_N ... + FAR_M ... + Ns Ms lambdas sigmas +; + +B = 10e5; +T_p = 1e-5; +f_s = B; +T_s = 1 / f_s; +f_c = 1e9; +duty_cycle_inv = 10; +T_r = T_p * duty_cycle_inv; +PRF = 1 / T_r; + +t_single_pulse = 0: 1 / f_s: T_r - 1 / f_s; + +c = 3e8; + +N = length(t_single_pulse); +N_high = T_p * f_s; +M = round(T_p * B); +R_0 = 0; + +delta_R_wide = c ./ 2 ./ B; % 宽带情况下的距离分辨力 +delta_R_narrow = T_p .* c ./ 2; % 窄带情况下的距离分辨力 + +N_wide = FAR_N; % 宽带情况下，发射 480 个脉冲 + +N_narrow = round(N_wide / M); + +ranges_wide = R_0 + (0:N_wide) * delta_R_wide; % [1000, 4000] +ranges_narrow = R_0 + (0:N_narrow) * delta_R_narrow; % [1000, 4000] + +% betas_wide = (linspace(1, 0.1, FAR_M * FAR_N) + 1j * linspace(0.1, 1, FAR_M * FAR_N))'; +betas_wide = zeros(FAR_M * FAR_N, 1) + 1; +betas_narrow = zeros(FAR_N, 1); +for i = 1: FAR_N + for j = 1: FAR_M + idx = (i-1) * FAR_M + j; + betas_narrow(i) = betas_narrow(i) + betas_wide(idx) * exp(1j * 2 * pi * 2 * delta_R_wide * j / c); + end +end + +trail_times = 250; +% method = "debiased_LASSO"; +% method = "LASSO"; +% method = "BP"; +method = "cVAMPro"; +% method = "debiased_LASSO_FISTA"; + +LASSO_lambda = 0.1; +tau = 1e-6; +iter_max = 1000; + +Ns = [256]; +Ms = [16]; +sigmas = [0.01]; +lambdas = zeros(length(Ns), 1); +for i = 1: length(Ns) + lambdas(i) = Map(Ns(i) + "_" + Ms(i)); +end diff --git a/wide vs narrow/single_point_test/query2.m b/wide vs narrow/single_point_test/query2.m new file mode 100755 index 0000000..403d23e --- /dev/null +++ b/wide vs narrow/single_point_test/query2.m @@ -0,0 +1,39 @@ +function [recovery_results, sigma_w2, thresholds, C_n, A, Lambda, Lambda_C, x] = query2(FAR_N, FAR_M, sigma, LASSO_lambda) + global trail_times method; + + if FAR_N == 2 + lambda_length = 1; + else + lambda_length = FAR_M; + end + + filename = ... + "./data3/" + ... + "FAR" + "_" + ... + string(FAR_N) + "_" + ... + string(FAR_M)+ "_" + ... + method + "_" + ... + string(sigma) + "_" + ... + string(lambda_length) + "_" + ... + string(LASSO_lambda) + "_" + ... + trail_times + ... + ".mat"; + + filename = ... + "./data3/" + ... + "FAR" + "_" + ... + string(FAR_N) + "_" + ... + string(FAR_M)+ "_" + ... + method + "_" + ... + string(sigma) + "_" + ... + trail_times + ... + ".mat"; + + load( ... + filename, ... + "recovery_results", "sigma_w2", "thresholds",... + "C_n", "A", "Lambda", "Lambda_C", "x" ... + ); + + return; +end \ No newline at end of file diff --git a/wide vs narrow/single_point_test/recovery.m b/wide vs narrow/single_point_test/recovery.m new file mode 100755 index 0000000..faa2b04 --- /dev/null +++ b/wide vs narrow/single_point_test/recovery.m @@ -0,0 +1,52 @@ +function x_hat = recovery(A, y_noise, method) + if method == "debiased_LASSO" + sz = size(A); + N = sz(2); + LASSO_lambda = 0.1; + + gamma = sz(1) / sz(2); + + cvx_begin quiet + variable x_LASSO(N) complex + minimize(LASSO_lambda * norm(x_LASSO, 1) + norm(y_noise - A * x_LASSO, 2)) + cvx_end + + rho_active = sum(abs(x_LASSO) > 1e-3)/N; + Q_hat = (gamma - rho_active)/(1 - rho_active); + Rho = sum((abs(x_LASSO) > 1e-3).* (2 - LASSO_lambda./(Q_hat*abs(x_LASSO) + LASSO_lambda))) / 2 / N; + diff = 1; + while(diff > 1e-4) + Rho_pre = Rho; + Rho = sum((abs(x_LASSO) > 1e-3).* (2 - LASSO_lambda./((gamma-Rho)/(1-Rho)*abs(x_LASSO) + LASSO_lambda))) / 2 / N; + diff = abs(Rho - Rho_pre); + end + Q_hat = (gamma-Rho)/(1-Rho); + x_d_CROD = x_LASSO + A'*(y_noise - A*x_LASSO)/Q_hat; + x_hat = x_d_CROD; + + elseif method == "LASSO" + sz = size(A); + N = sz(2); + LASSO_lambda = 0.1; + cvx_begin quiet + variable x_LASSO(N) complex + minimize(LASSO_lambda * norm(x_LASSO, 1) + norm(y_noise - A * x_LASSO, 2)) + cvx_end + + x_hat = x_LASSO; + + elseif method == "BP" + sz = size(A); + N = sz(2); + cvx_begin quiet + variable x_hat(N) complex + minimize(norm(x_hat, 1)) + subject to + A * x_hat == y_noise + cvx_end + else + global lambda tau iter_max; + [x_hat_wl, x_hat] = cVAMPro(y_noise, A, lambda, tau, iter_max); + end +end + diff --git a/wide vs narrow/single_point_test/signle_point_mento_carlo_test.m b/wide vs narrow/single_point_test/signle_point_mento_carlo_test.m new file mode 100755 index 0000000..e848d4f --- /dev/null +++ b/wide vs narrow/single_point_test/signle_point_mento_carlo_test.m @@ -0,0 +1,269 @@ +% 单目标检测点测试 + +clc; clear; close all; + +parameters; + +% Ns = [4, 8, 16, 32, 64, 128]; + +TEST_1 = true; +TEST_2 = false; +TEST_3 = false; +TEST_4 = false; + +if TEST_1 + for N_idx = 1: length(Ns) + for M_idx = 1: length(Ms) + for sigma_idx = 1: length(sigmas) + FAR_N = Ns(N_idx); + FAR_M = Ms(M_idx); + sigma = sigmas(sigma_idx); + for lambda_idx = 1: length(lambdas) + LASSO_lambda = lambdas(lambda_idx); + + % epi = 0; + % [C_n, A] = get_Psi(FAR_N, FAR_M, epi); + % fac = A * A'; + % A = A ./ sqrt(abs(fac(1, 1))); + % f_n = f_c + C_n * B / M; + % + % [Lambda, Lambda_C, x] = get_sparse_vector(FAR_N, FAR_M, betas_wide, true); + % [recovery_results, sigma_w2, thresholds] = All_Recovery(A, x, LASSO_lambda, sigma); + + [recovery_results, sigma_w2, thresholds, C_n, A, Lambda, Lambda_C, x] = query2(FAR_N, FAR_M, sigma, LASSO_lambda); + + i1 = 1; i2 = 1; + Lambda_distributes = zeros(2, length(Lambda) * trail_times); + Lambda_C_distributes = zeros(2, length(Lambda_C) * trail_times); + x_1_add = length(length(Lambda) + length(Lambda_C)); + + for t = 1:trail_times + x_0_hat = squeeze(recovery_results(1, t, :)); + x_1_hat = squeeze(recovery_results(2, t, :)) - x; + if anynan(x_1_hat) + continue; + end + + for i = 1: length(Lambda) + length(Lambda_C) + if ismember(i, Lambda_C) + Lambda_C_distributes(1, i1) = x_0_hat(i); + Lambda_C_distributes(2, i1) = x_1_hat(i); + i1 = i1 + 1; + elseif ismember(i, Lambda) + Lambda_distributes(1, i2) = x_0_hat(i); + Lambda_distributes(2, i2) = x_1_hat(i); + i2 = i2 + 1; + end + end + end + real_H00 = real(Lambda_C_distributes(1, 1:i1-1)); + real_H01 = real(Lambda_distributes(1, 1:i2-1)); + real_H10 = real(Lambda_C_distributes(2, 1:i1-1)); + real_H11 = real(Lambda_distributes(2, 1:i2-1)); + + figure; + subplot(2, 2, 1); histfit(real_H00); xlabel("x\_hat"); ylabel("times"); + subplot(2, 2, 2); histfit(real_H01); xlabel("x\_hat"); ylabel("times"); + subplot(2, 2, 3); histfit(real_H10); xlabel("x\_hat"); ylabel("times"); + subplot(2, 2, 4); histfit(real_H11); xlabel("x\_hat"); ylabel("times"); + + % figure; + % data = real(Lambda_C_distributes(2, :)); + % xx = linspace(min(data), max(data), 1e3); + % yy = normpdf(x, mean(data), std(data)); + % hold on; + % histfit(data); + + xlabel("x\_hat"); + ylabel("times"); + legend("lambda = " + string(LASSO_lambda) + ", sigma = " + string(sigma)); + + means = [ ... + mean(Lambda_C_distributes(1, 1:i1-1)), ... + mean(Lambda_distributes(1, i2-1)), ... + mean(Lambda_C_distributes(2, 1:i1-1)), ... + mean(Lambda_distributes(2, i2-1)), ... + mean([Lambda_distributes(1, i2-1) Lambda_C_distributes(1, 1:i1-1)]) ... + ]; + + stds = [ ... + std(Lambda_C_distributes(1, 1:i1-1)), ... + std(Lambda_distributes(1, i2-1)), ... + std(Lambda_C_distributes(2, 1:i1-1)), ... + std(Lambda_distributes(2, i2-1)), ... + std([Lambda_distributes(1, i2-1) Lambda_C_distributes(1, 1:i1-1)]) ... + ]; + + + fprintf("H_00: mu = %.15f, std = %.15f\n", means(1), stds(1)); + fprintf("H_01: mu = %.15f, std = %.15f\n", means(2), stds(2)); + fprintf("H_10: mu = %.15f, std = %.15f\n", means(3), stds(3)); + fprintf("H_11: mu = %.15f, std = %.15f\n", means(4), stds(4)); + fprintf("H_0: mu = %.15f, std = %.15f\n", means(5), stds(5)); + + [ss, vs, H0_mean, H0_max, H1_mean, H1_max] = get_sigma_var(recovery_results, sigma_w2, x); + end + end + end + end +end + +if TEST_2 + lg_lambdas = -5: 0.01: -0.75; + % lg_lambdas = -3: 0.01: -2.5; + lambdas = 10 .^ lg_lambdas; + LASSO_x_hat_maps = zeros(length(Ns), length(Ms), length(sigmas), length(lambdas)); + LASSO_lambdas = zeros(length(Ns), length(Ms), length(sigmas)); + for sigma_idx = 1: length(c) + sigma = sigmas(sigma_idx); + figure; + for N_idx = 1: length(Ns) + FAR_N = Ns(N_idx); + for M_idx = 1: length(Ms) + FAR_M = Ms(M_idx); + subplot(length(Ns), length(Ms), (N_idx-1)*length(Ms)+M_idx); + epi = 0; + [C_n, A] = get_Psi(FAR_N, FAR_M, epi); + fac = A * A'; + A = A ./ sqrt(abs(fac(1, 1))); + f_n = f_c + C_n * B / M; + + [Lambda, Lambda_C, x] = get_sparse_vector(FAR_N, FAR_M, betas_wide, true); + [LASSO_lambda, LASSO_x_hat_map] = get_lambda(A, x, FAR_N, FAR_M, sigma, lambdas); + LASSO_x_hat_maps(N_idx, M_idx, sigma_idx, :) = LASSO_x_hat_map; + LASSO_lambdas(N_idx, M_idx, sigma_idx) = LASSO_lambda; + + semilogx(lambdas, LASSO_x_hat_map / FAR_M); + xlabel("\lambda in LASSO"); + ylabel("mean(x\_hat(\Lambda))"); + yline(1); + ylim([-0.2, 1.2]); + legend("N = " + string(FAR_N), "M = " + string(FAR_M)); + end + end + end + save("tmp.mat", lambdas, LASSO_x_hat_maps, LASSO_lambdas); +end + +if TEST_3 + for N_idx = 1: length(Ns) + for M_idx = 1: length(Ms) + for sigma_idx = 1: length(sigmas) + FAR_N = Ns(N_idx); + FAR_M = Ms(M_idx); + sigma = sigmas(sigma_idx); + for lambda_idx = 1: length(lambdas) + LASSO_lambda = lambdas(lambda_idx); + + % epi = 0; + % [C_n, A] = get_Psi(FAR_N, FAR_M, epi); + % fac = A * A'; + % A = A ./ sqrt(abs(fac(1, 1))); + % f_n = f_c + C_n * B / M; + % + % [Lambda, Lambda_C, x] = get_sparse_vector(FAR_N, FAR_M, betas_wide, true); + % [recovery_results, sigma_w2, thresholds] = All_Recovery(A, x, LASSO_lambda, sigma); + + [recovery_results, sigma_w2, thresholds, C_n, A, Lambda, Lambda_C, x] = query2(FAR_N, FAR_M, sigma, LASSO_lambda); + + i1 = 1; i2 = 1; + Lambda_distributes = zeros(2, length(Lambda) * trail_times); + Lambda_C_distributes = zeros(2, length(Lambda_C) * trail_times); + x_1_add = length(length(Lambda) + length(Lambda_C)); + + for t = 1:trail_times + x_0_hat = squeeze(recovery_results(1, t, :)); + x_1_hat = squeeze(recovery_results(2, t, :)) - x; + + for i = 1: length(Lambda) + length(Lambda_C) + if ismember(i, Lambda_C) + Lambda_C_distributes(1, i1) = x_0_hat(i); + Lambda_C_distributes(2, i1) = x_1_hat(i); + i1 = i1 + 1; + elseif ismember(i, Lambda) + Lambda_distributes(1, i2) = x_0_hat(i); + Lambda_distributes(2, i2) = x_1_hat(i); + i2 = i2 + 1; + end + end + end + real_H00 = real(Lambda_C_distributes(1, 1:i1-1)); + real_H01 = real(Lambda_distributes(1, 1:i2-1)); + real_H10 = real(Lambda_C_distributes(2, 1:i1-1)); + real_H11 = real(Lambda_distributes(2, 1:i2-1)); + + + figure; + subplot(2, 2, 1); histfit(real_H00); xlabel("x\_hat"); ylabel("times"); + subplot(2, 2, 2); histfit(real_H01); xlabel("x\_hat"); ylabel("times"); + subplot(2, 2, 3); histfit(real_H10); xlabel("x\_hat"); ylabel("times"); + subplot(2, 2, 4); histfit(real_H11); xlabel("x\_hat"); ylabel("times"); + + figure; + data = real(Lambda_C_distributes(2, :)); + xx = linspace(min(data), max(data), 1e3); + yy = normpdf(x, mean(data), std(data)); + hold on; + histfit(data); + + xlabel("x\_hat"); + ylabel("times"); + legend("lambda = " + string(LASSO_lambda) + ", sigma = " + string(sigma)); + + means = [ ... + mean(Lambda_C_distributes(1, :)), ... + mean(Lambda_distributes(1, :)), ... + mean(Lambda_C_distributes(2, :)), ... + mean(Lambda_distributes(2, :)), ... + mean([Lambda_distributes(1, :) Lambda_C_distributes(1, :)]), ... + mean([Lambda_distributes(1, :) Lambda_C_distributes(1, :) Lambda_C_distributes(2, :)]), ... + mean([Lambda_distributes(2, :)]) ... + ]; + + stds = [ ... + std(Lambda_C_distributes(1, :)), ... + std(Lambda_distributes(1, :)), ... + std(Lambda_C_distributes(2, :)), ... + std(Lambda_distributes(2, :)), ... + std([Lambda_distributes(1, :) Lambda_C_distributes(1, :)]), ... + std([Lambda_distributes(1, :) Lambda_C_distributes(1, :) Lambda_C_distributes(2, :)]), ... + std([Lambda_distributes(2, :)]) ... + ]; + + + fprintf("H_00: mu = %.15f, std = %.15f\n", means(1), stds(1)); + fprintf("H_01: mu = %.15f, std = %.15f\n", means(2), stds(2)); + fprintf("H_10: mu = %.15f, std = %.15f\n", means(3), stds(3)); + fprintf("H_11: mu = %.15f, std = %.15f\n", means(4), stds(4)); + fprintf("H_0: mu = %.15f, std = %.15f\n", means(5), stds(5)); + fprintf("H_0 + Lambda^C: mu = %.15f, std = %.15f\n", means(6), stds(6)); + fprintf("H_1 Lambda: mu = %.15f, std = %.15f\n\n", means(7), stds(7)); + + [ss, vs, H0_mean, H0_max, H1_mean, H1_max] = get_sigma_var(recovery_results, sigma_w2, x); + end + end + end + end +end + + +% for N_idx = 1: length(Ns) +% for M_idx = 1: length(Ms) +% for lambda_idx = 1: length(lambdas) +% for sigma_idx = 1: length(sigmas) +% FAR_N = Ns(N_idx); +% FAR_M = Ms(M_idx); +% LASSO_lambda = lambdas(lambda_idx); +% sigma = sigmas(sigma_idx); +% [recovery_results, sigma_w2, thresholds, C_n, A, Lambda, Lambda_C, x] = query(FAR_N, FAR_M, sigma, LASSO_lambda); +% end +% end +% end +% end + + +% figure; +% subplot(2, 2, 1); histfit(real(Lambda_C_distributes(1, :))); xlabel("x\_hat"); ylabel("times"); +% subplot(2, 2, 2); histfit(real(Lambda_distributes(1, :))); xlabel("x\_hat"); ylabel("times"); +% subplot(2, 2, 3); histfit(real(Lambda_C_distributes(2, :))); xlabel("x\_hat"); ylabel("times"); +% subplot(2, 2, 4); histfit(real(Lambda_distributes(2, :))); xlabel("x\_hat"); ylabel("times"); \ No newline at end of file diff --git a/wide vs narrow/single_point_test/signle_point_mento_carlo_train.m b/wide vs narrow/single_point_test/signle_point_mento_carlo_train.m new file mode 100755 index 0000000..aab72f0 --- /dev/null +++ b/wide vs narrow/single_point_test/signle_point_mento_carlo_train.m @@ -0,0 +1,53 @@ +% 单目标检测点测试 + +clc; clear; close all; + +parameters; + +for N_idx = 1: length(Ns) + for M_idx = 1: length(Ms) + for sigma_idx = 1: length(sigmas) + FAR_N = Ns(N_idx); + FAR_M = Ms(M_idx); + % LASSO_lambda = lambdas(lambda_idx); + sigma = sigmas(sigma_idx); + for lambda_idx = 1: length(lambdas) + LASSO_lambda = lambdas(lambda_idx); + + epi = 0; + Psi_filename = "./results/Signal_Model_" + string(FAR_N) + "_" + string(FAR_M) + ".mat"; + load(Psi_filename, "C_n", "A", "ref_lambda"); + LASSO_lambda = ref_lambda; + % [C_n, A] = get_Psi(FAR_N, FAR_M, epi); + f_n = f_c + C_n * B / M; + + [Lambda, Lambda_C, x] = get_sparse_vector(FAR_N, FAR_M, betas_wide, true); + filename = ... + "./data3/" + ... + "FAR" + "_" + ... + string(FAR_N) + "_" + ... + string(FAR_M)+ "_" + ... + method + "_" + ... + string(sigma) + "_" + ... + trail_times + ... + ".mat"; + if exist(filename, "file") + load( ... + filename, ... + "recovery_results", "sigma_w2", "thresholds",... + "C_n", "A", "Lambda", "Lambda_C", "x" ... + ); + else + % Mento Carlo Recovery + [recovery_results, sigma_w2, thresholds] = All_Recovery(A, x, LASSO_lambda, sigma); + + save( ... + filename, ... + "recovery_results", "sigma_w2", "thresholds",... + "C_n", "A", "Lambda", "Lambda_C", "x" ... + ); + end + end + end + end +end \ No newline at end of file diff --git a/wide vs narrow/single_point_test/simulator.m b/wide vs narrow/single_point_test/simulator.m new file mode 100644 index 0000000..ebaa973 --- /dev/null +++ b/wide vs narrow/single_point_test/simulator.m @@ -0,0 +1,196 @@ +%% Initial +clc; clear; close all; + +trail_times = 250; +method = "cVAMPro"; +Ns = [64, 128, 256, 512, 1024]; +Ms = [4, 8, 16, 32]; +sigma = 0.01; +tau = 1e-6; +iter_max = 100; +epi = 0; + +lg_lambdas = -5: 0.05: -1; +lambdas = 10 .^ lg_lambdas; + +H0_REE_means = zeros(length(Ns), length(Ms)); +H0_REE_maxs = zeros(length(Ns), length(Ms)); +H1_REE_means = zeros(length(Ns), length(Ms)); +H1_REE_maxs = zeros(length(Ns), length(Ms)); + +figure; + +for N_idx = 1: length(Ns) + for M_idx = 1: length(Ms) + subplot(length(Ns), length(Ms), (N_idx - 1) * length(Ms) + M_idx); + FAR_N = Ns(N_idx); + FAR_M = Ms(M_idx); + + fprintf("Simluating: N = %d, M = %d\n\n", FAR_N, FAR_M); + + MN = FAR_N * FAR_M; + betas_wide = zeros(FAR_M * FAR_N, 1) + 1; + [Lambda, Lambda_C, x] = get_sparse_vector(FAR_N, FAR_M, betas_wide, true); + + % get lambda + signal_model_filename = "./Signal_Model/Signal_Model_" + string(FAR_N) + "_" + string(FAR_M) + ".mat"; + if exist(signal_model_filename, "file") + load(signal_model_filename, "A", "C_n", "ref_lambda"); + else + curr_metric = FAR_N; + for TT = 1: 20 + [C_n, A] = get_Psi(FAR_N, FAR_M, epi); + + noise = get_noise(sigma, FAR_N, 1); + y_noise = A * x + noise; + + MSEs = zeros(length(lambdas), 1); + all_x_hat = zeros(length(x), length(lambdas)); + + for lambda_idx = 1: length(lambdas) + lambda = lambdas(lambda_idx); + [x_LASSO, x_hat_d] = cVAMPro(y_noise, A, lambda, tau, iter_max); + MSEs(lambda_idx) = sum(real(x_LASSO - x)); + all_x_hat(:, lambda_idx) = x_LASSO; + end + + [metric, idx] = min(abs(MSEs)); + if curr_metric > metric + ref_lambda = lambdas(idx); + save(signal_model_filename, "A", "C_n", "ref_lambda"); + curr_metric = metric; + end + end + end + fprintf("get lambda: lambda = %.15f. \n", ref_lambda); + fprintf("Signal Model filename: " + signal_model_filename + "\n\n"); + + % train + LASSO_lambda = ref_lambda; + simulate_results_filename = ... + "./data5/" + ... + "FAR" + "_" + ... + string(FAR_N) + "_" + ... + string(FAR_M)+ "_" + ... + method + "_" + ... + string(sigma) + "_" + ... + trail_times + ... + ".mat"; + if exist(simulate_results_filename, "file") + load( ... + simulate_results_filename, ... + "recovery_results", "sigma_w2", "thresholds",... + "C_n", "A", "Lambda", "Lambda_C", "x" ... + ); + else + % Mento Carlo Recovery + [recovery_results, sigma_w2, thresholds] = All_Recovery2(A, x, LASSO_lambda, sigma, trail_times, LASSO_lambda, tau); + + save( ... + simulate_results_filename, ... + "recovery_results", "sigma_w2", "thresholds",... + "C_n", "A", "Lambda", "Lambda_C", "x" ... + ); + end + fprintf("Simulate results filename: " + simulate_results_filename + "\n\n"); + + % test + i1 = 1; i2 = 1; + Lambda_distributes = zeros(2, length(Lambda) * trail_times); + Lambda_C_distributes = zeros(2, length(Lambda_C) * trail_times); + + for t = 1:trail_times + x_0_hat = squeeze(recovery_results(1, t, :)); + x_1_hat = squeeze(recovery_results(2, t, :)) - x; + if anynan(x_1_hat) + continue; + end + + for i = 1: length(Lambda) + length(Lambda_C) + if ismember(i, Lambda_C) + Lambda_C_distributes(1, i1) = x_0_hat(i); + Lambda_C_distributes(2, i1) = x_1_hat(i); + i1 = i1 + 1; + elseif ismember(i, Lambda) + Lambda_distributes(1, i2) = x_0_hat(i); + Lambda_distributes(2, i2) = x_1_hat(i); + i2 = i2 + 1; + end + end + end + real_H00 = real(Lambda_C_distributes(1, 1:i1-1)); + real_H01 = real(Lambda_distributes(1, 1:i2-1)); + real_H10 = real(Lambda_C_distributes(2, 1:i1-1)); + real_H11 = real(Lambda_distributes(2, 1:i2-1)); + + % figure; + % subplot(2, 2, 1); histfit(real_H00); xlabel("x\_hat"); ylabel("times"); title("H_0 (not in support set)") + % subplot(2, 2, 2); histfit(real_H01); xlabel("x\_hat"); ylabel("times"); title("H_0 (in support set)") + % subplot(2, 2, 3); histfit(real_H10); xlabel("x\_hat"); ylabel("times"); title("H_1 (not in support set)") + % subplot(2, 2, 4); histfit(real_H11); xlabel("x\_hat"); ylabel("times"); title("H_1 (in support set)") + % sgtitle("N = " + string(FAR_N) + ", M = " + string(FAR_M)); + + histfit(real_H11); xlabel("x\_hat"); ylabel("times"); title("N = " + string(FAR_N) + ", M = " + string(FAR_M)); + + means = [ ... + mean(Lambda_C_distributes(1, 1:i1-1)), ... + mean(Lambda_distributes(1, i2-1)), ... + mean(Lambda_C_distributes(2, 1:i1-1)), ... + mean(Lambda_distributes(2, i2-1)), ... + mean([Lambda_distributes(1, i2-1) Lambda_C_distributes(1, 1:i1-1)]) ... + ]; + + stds = [ ... + std(Lambda_C_distributes(1, 1:i1-1)), ... + std(Lambda_distributes(1, i2-1)), ... + std(Lambda_C_distributes(2, 1:i1-1)), ... + std(Lambda_distributes(2, i2-1)), ... + std([Lambda_distributes(1, i2-1) Lambda_C_distributes(1, 1:i1-1)]) ... + ]; + + + fprintf("H_00: mu = %.15f, std = %.15f\n", means(1), stds(1)); + fprintf("H_01: mu = %.15f, std = %.15f\n", means(2), stds(2)); + fprintf("H_10: mu = %.15f, std = %.15f\n", means(3), stds(3)); + fprintf("H_11: mu = %.15f, std = %.15f\n", means(4), stds(4)); + fprintf("H_0: mu = %.15f, std = %.15f\n", means(5), stds(5)); + + [ss, vs, H0_mean, H0_max, H1_mean, H1_max] = get_sigma_var(recovery_results, sigma_w2, x); + H0_REE_means(N_idx, M_idx) = H0_mean; + H0_REE_maxs(N_idx, M_idx) = H0_max; + H1_REE_means(N_idx, M_idx) = H1_mean; + H1_REE_maxs(N_idx, M_idx) = H1_max; + fprintf("Test Complete \n\n"); + end +end +sgtitle("H_1 (in support set)"); + +figure; + +subplot(221); +h = heatmap(Ns, Ms, H0_REE_means'); +h.XLabel = "N"; +h.YLabel = "M"; +h.Title = "H_0, mean(REE)"; + +subplot(222); +h = heatmap(Ns, Ms, H0_REE_maxs'); +h.XLabel = "N"; +h.YLabel = "M"; +h.Title = "H_0, max(REE)"; + +subplot(223); +h = heatmap(Ns, Ms, H1_REE_means'); +h.XLabel = "N"; +h.YLabel = "M"; +h.Title = "H_1, mean(REE)"; + +subplot(224); +h = heatmap(Ns, Ms, H1_REE_maxs'); +h.XLabel = "N"; +h.YLabel = "M"; +h.Title = "H_1, max(REE)"; + + + + diff --git a/wide vs narrow/single_point_test/swtest.m b/wide vs narrow/single_point_test/swtest.m new file mode 100755 index 0000000..5112e5a --- /dev/null +++ b/wide vs narrow/single_point_test/swtest.m @@ -0,0 +1,273 @@ +function [H, pValue, W] = swtest(x, alpha) +%SWTEST Shapiro-Wilk parametric hypothesis test of composite normality. +% [H, pValue, SWstatistic] = SWTEST(X, ALPHA) performs the +% Shapiro-Wilk test to determine if the null hypothesis of +% composite normality is a reasonable assumption regarding the +% population distribution of a random sample X. The desired significance +% level, ALPHA, is an optional scalar input (default = 0.05). +% +% The Shapiro-Wilk and Shapiro-Francia null hypothesis is: +% "X is normal with unspecified mean and variance." +% +% This is an omnibus test, and is generally considered relatively +% powerful against a variety of alternatives. +% Shapiro-Wilk test is better than the Shapiro-Francia test for +% Platykurtic sample. Conversely, Shapiro-Francia test is better than the +% Shapiro-Wilk test for Leptokurtic samples. +% +% When the series 'X' is Leptokurtic, SWTEST performs the Shapiro-Francia +% test, else (series 'X' is Platykurtic) SWTEST performs the +% Shapiro-Wilk test. +% +% [H, pValue, SWstatistic] = SWTEST(X, ALPHA) +% +% Inputs: +% X - a vector of deviates from an unknown distribution. The observation +% number must exceed 3 and less than 5000. +% +% Optional inputs: +% ALPHA - The significance level for the test (default = 0.05). +% +% Outputs: +% SWstatistic - The test statistic (non normalized). +% +% pValue - is the p-value, or the probability of observing the given +% result by chance given that the null hypothesis is true. Small values +% of pValue cast doubt on the validity of the null hypothesis. +% +% H = 0 => Do not reject the null hypothesis at significance level ALPHA. +% H = 1 => Reject the null hypothesis at significance level ALPHA. +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Copyright (c) 17 March 2009 by Ahmed Ben Sada % +% Department of Finance, IHEC Sousse - Tunisia % +% Email: ahmedbensaida@yahoo.com % +% $ Revision 3.0 $ Date: 18 Juin 2014 $ % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + +% +% References: +% +% - Royston P. "Remark AS R94", Applied Statistics (1995), Vol. 44, +% No. 4, pp. 547-551. +% AS R94 -- calculates Shapiro-Wilk normality test and P-value +% for sample sizes 3 <= n <= 5000. Handles censored or uncensored data. +% Corrects AS 181, which was found to be inaccurate for n > 50. +% Subroutine can be found at: http://lib.stat.cmu.edu/apstat/R94 +% +% - Royston P. "A pocket-calculator algorithm for the Shapiro-Francia test +% for non-normality: An application to medicine", Statistics in Medecine +% (1993a), Vol. 12, pp. 181-184. +% +% - Royston P. "A Toolkit for Testing Non-Normality in Complete and +% Censored Samples", Journal of the Royal Statistical Society Series D +% (1993b), Vol. 42, No. 1, pp. 37-43. +% +% - Royston P. "Approximating the Shapiro-Wilk W-test for non-normality", +% Statistics and Computing (1992), Vol. 2, pp. 117-119. +% +% - Royston P. "An Extension of Shapiro and Wilk's W Test for Normality +% to Large Samples", Journal of the Royal Statistical Society Series C +% (1982a), Vol. 31, No. 2, pp. 115-124. +% + +% +% Ensure the sample data is a VECTOR. +% + +if numel(x) == length(x) + x = x(:); % Ensure a column vector. +else + error(' Input sample ''X'' must be a vector.'); +end + +% +% Remove missing observations indicated by NaN's and check sample size. +% + +x = x(~isnan(x)); + +if length(x) < 3 + error(' Sample vector ''X'' must have at least 3 valid observations.'); +end + +if length(x) > 5000 + warning('Shapiro-Wilk test might be inaccurate due to large sample size ( > 5000).'); +end + +if max(x) == min(x) + H = 0; + pValue = 0; + W = 0; + return; +end + +% +% Ensure the significance level, ALPHA, is a +% scalar, and set default if necessary. +% + +if (nargin >= 2) && ~isempty(alpha) + if ~isscalar(alpha) + error(' Significance level ''Alpha'' must be a scalar.'); + end + if (alpha <= 0 || alpha >= 1) + error(' Significance level ''Alpha'' must be between 0 and 1.'); + end +else + alpha = 0.05; +end + +% First, calculate the a's for weights as a function of the m's +% See Royston (1992, p. 117) and Royston (1993b, p. 38) for details +% in the approximation. + +x = sort(x); % Sort the vector X in ascending order. +n = length(x); +mtilde = norminv(((1:n)' - 3/8) / (n + 1/4)); +weights = zeros(n,1); % Preallocate the weights. + +if kurtosis(x) > 3 + + % The Shapiro-Francia test is better for leptokurtic samples. + + weights = 1/sqrt(mtilde'*mtilde) * mtilde; + + % + % The Shapiro-Francia statistic W' is calculated to avoid excessive + % rounding errors for W' close to 1 (a potential problem in very + % large samples). + % + + W = (weights' * x)^2 / ((x - mean(x))' * (x - mean(x))); + + % Royston (1993a, p. 183): + nu = log(n); + u1 = log(nu) - nu; + u2 = log(nu) + 2/nu; + mu = -1.2725 + (1.0521 * u1); + sigma = 1.0308 - (0.26758 * u2); + + newSFstatistic = log(1 - W); + + % + % Compute the normalized Shapiro-Francia statistic and its p-value. + % + + NormalSFstatistic = (newSFstatistic - mu) / sigma; + + % Computes the p-value, Royston (1993a, p. 183). + pValue = 1 - normcdf(real(NormalSFstatistic), 0, 1); + +else + + % The Shapiro-Wilk test is better for platykurtic samples. + + c = 1/sqrt(mtilde'*mtilde) * mtilde; + u = 1/sqrt(n); + + % Royston (1992, p. 117) and Royston (1993b, p. 38): + PolyCoef_1 = [-2.706056 , 4.434685 , -2.071190 , -0.147981 , 0.221157 , c(n)]; + PolyCoef_2 = [-3.582633 , 5.682633 , -1.752461 , -0.293762 , 0.042981 , c(n-1)]; + + % Royston (1992, p. 118) and Royston (1993b, p. 40, Table 1) + PolyCoef_3 = [-0.0006714 , 0.0250540 , -0.39978 , 0.54400]; + PolyCoef_4 = [-0.0020322 , 0.0627670 , -0.77857 , 1.38220]; + PolyCoef_5 = [0.00389150 , -0.083751 , -0.31082 , -1.5861]; + PolyCoef_6 = [0.00303020 , -0.082676 , -0.48030]; + + PolyCoef_7 = [0.459 , -2.273]; + + weights(n) = polyval(PolyCoef_1 , u); + weights(1) = -weights(n); + + if n > 5 + weights(n-1) = polyval(PolyCoef_2 , u); + weights(2) = -weights(n-1); + + count = 3; + phi = (mtilde'*mtilde - 2 * mtilde(n)^2 - 2 * mtilde(n-1)^2) / ... + (1 - 2 * weights(n)^2 - 2 * weights(n-1)^2); + else + count = 2; + phi = (mtilde'*mtilde - 2 * mtilde(n)^2) / ... + (1 - 2 * weights(n)^2); + end + + % Special attention when n = 3 (this is a special case). + if n == 3 + % Royston (1992, p. 117) + weights(1) = 1/sqrt(2); + weights(n) = -weights(1); + phi = 1; + end + + % + % The vector 'WEIGHTS' obtained next corresponds to the same coefficients + % listed by Shapiro-Wilk in their original test for small samples. + % + + weights(count : n-count+1) = mtilde(count : n-count+1) / sqrt(phi); + + % + % The Shapiro-Wilk statistic W is calculated to avoid excessive rounding + % errors for W close to 1 (a potential problem in very large samples). + % + + W = (weights' * x) ^2 / ((x - mean(x))' * (x - mean(x))); + + % + % Calculate the normalized W and its significance level (exact for + % n = 3). Royston (1992, p. 118) and Royston (1993b, p. 40, Table 1). + % + + newn = log(n); + + if (n >= 4) && (n <= 11) + + mu = polyval(PolyCoef_3 , n); + sigma = exp(polyval(PolyCoef_4 , n)); + gam = polyval(PolyCoef_7 , n); + + newSWstatistic = -log(gam-log(1-W)); + + elseif n > 11 + + mu = polyval(PolyCoef_5 , newn); + sigma = exp(polyval(PolyCoef_6 , newn)); + + newSWstatistic = log(1 - W); + + elseif n == 3 + mu = 0; + sigma = 1; + newSWstatistic = 0; + end + + % + % Compute the normalized Shapiro-Wilk statistic and its p-value. + % + + NormalSWstatistic = (newSWstatistic - mu) / sigma; + + % NormalSWstatistic is referred to the upper tail of N(0,1), + % Royston (1992, p. 119). + pValue = 1 - normcdf(NormalSWstatistic, 0, 1); + + % Special attention when n = 3 (this is a special case). + if n == 3 + pValue = 6/pi * (asin(sqrt(W)) - asin(sqrt(3/4))); + % Royston (1982a, p. 121) + end + +end + +% +% To maintain consistency with existing Statistics Toolbox hypothesis +% tests, returning 'H = 0' implies that we 'Do not reject the null +% hypothesis at the significance level of alpha' and 'H = 1' implies +% that we 'Reject the null hypothesis at significance level of alpha.' +% + +H = (alpha >= pValue); \ No newline at end of file diff --git a/wide vs narrow/single_point_test/test_lambda.m b/wide vs narrow/single_point_test/test_lambda.m new file mode 100755 index 0000000..d2516af --- /dev/null +++ b/wide vs narrow/single_point_test/test_lambda.m @@ -0,0 +1,102 @@ +clc; clear; close all; + +N = 256; +M = 16; +MN = N * M; +sigmas = [0.01]; +% noise_sigma = 0.01; + +filename = "./results/Signal_Model_" + string(N) + "_" + string(M) + ".mat"; +curr_metric = N; +for TT = 1: 100 + TT +tau = 1e-6; +iter_max = 25; +epi = 0; +[C_n, A] = get_Psi(N, M, epi); +betas_wide = zeros(M * N, 1) + 1; +[Lambda, Lambda_C, x] = get_sparse_vector(N, M, betas_wide, true); + +lg_lambdas = -3: 0.05: -1; +lambdas = 10 .^ lg_lambdas; +method = "cVAMPro"; + + +for sigma_idx = 1: length(sigmas) + sigma = sigmas(sigma_idx); + noise = get_noise(sigma, N, 1); + y_noise = A * x + noise; + + MSEs = zeros(length(lambdas), 1); + all_x_hat = zeros(length(x), length(lambdas)); + for lambda_idx = 1: length(lambdas) + lambda = lambdas(lambda_idx); + if method == "cVAMPro" + [x_LASSO, x_hat_d] = cVAMPro(y_noise, A, lambda, tau, iter_max); + elseif method == "debiased\_LASSO" + [x_LASSO, sigma_w_2, threshold] = debiased_LASSO(A, y_noise, 0, sigma^2, lambda); + elseif method == "LASSO" + cvx_begin quiet + variable x_LASSO(MN) complex + minimize(lambda * norm(x_LASSO, 1) + norm(y_noise - A * x_LASSO, 2)) + cvx_end + elseif method == "FISTA" + x_LASSO = FISTA(y_noise, A, LASSO_lambda, 1e-5); + elseif method == "debiased\_LASSO\_FISTA" + [x_LASSO, sigma_w_2, threshold] = debiased_LASSO_FISTA(A, y_noise, 0, noise_sigma^2, lambda); + elseif method == "BP" + sz = size(A); + N = sz(2); + cvx_begin quiet + variable x_LASSO(N) complex + minimize(norm(x_LASSO, 1)) + subject to + A * x_LASSO == y_noise + cvx_end + end + MSEs(lambda_idx) = sum(real(x_LASSO - x)); + all_x_hat(:, lambda_idx) = x_LASSO; + end + abs_MSEs = abs(MSEs); + [metric, idx] = min(abs_MSEs); + if curr_metric > metric + ref_lambda = lambdas(idx); + save(filename, "A", "C_n", "ref_lambda"); + curr_metric = metric; + end + if curr_metric < 0.1 + figure; + subplot(211); + semilogx(lambdas, MSEs); + yline(0); + ylim([-length(Lambda) - 0.1, 0.1]); + xlabel("\lambda"); + ylabel("MSE (x\_hat - x)"); + title("N = " + string(N) + ", M = " + string(M) + ", \sigma = " + string(sigma) + ", " + method); + subplot(212); + plot(C_n); + xlabel("Frequency code"); + ylabel("Times"); + break + end + % for i = 1: length(C_n) + % fprintf("%d, ", C_n(i)) + % end + % fprintf("%.15f\n", sum(MSEs)); + + % filename = ... + % string(N) + "_" + ... + % string(M) + "_" + ... + % string(sigma) + "_" + ... + % method + ".mat"; + % save(filename, "lambdas", "MSEs", "x", "all_x_hat"); + + % figure; + % semilogx(lambdas, xx); + % yline(0); + % ylim([-length(Lambda) - 0.1, 0.1]); + % xlabel("\lambda"); + % ylabel("MSE (x\_hat - x)"); + % title("N = " + string(FAR_N) + ", M = " + string(FAR_M) + ", \sigma = " + string(sigma) + ", " + method); +end +end