Ksyer/FAR_CS
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Update wide vs narrow code

Browse Source
This commit is contained in:
Ksyer
2024-07-22 21:17:57 +08:00
parent 13f660826d
commit 9cbe8e31ea
14 changed files with 1642 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files
wide vs narrow/Main3.m
Executable
+253
View File
@@ -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");
wide vs narrow/single_point_test/All_Recovery.m
Executable
+104
View File
@@ -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
wide vs narrow/single_point_test/All_Recovery2.m
Executable
+95
View File
@@ -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
wide vs narrow/single_point_test/get_lambda.m
Executable
+19
View File
@@ -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
wide vs narrow/single_point_test/get_sigma_var.m
Executable
+71
View File
@@ -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
wide vs narrow/single_point_test/get_sparse_vector.m
Executable
+28
View File
@@ -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
wide vs narrow/single_point_test/parameters.m
Executable
+88
View File
@@ -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
wide vs narrow/single_point_test/query2.m
Executable
+39
View File
@@ -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
wide vs narrow/single_point_test/recovery.m
Executable
+52
View File
@@ -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
wide vs narrow/single_point_test/signle_point_mento_carlo_test.m
Executable
+269
View File
@@ -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");
wide vs narrow/single_point_test/signle_point_mento_carlo_train.m
Executable
+53
View File
@@ -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
wide vs narrow/single_point_test/simulator.m
+196
View File
@@ -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)";
wide vs narrow/single_point_test/swtest.m
Executable
+273
View File
@@ -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);
wide vs narrow/single_point_test/test_lambda.m
Executable
+102
View File
@@ -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