function threshold = FAR_simu(P_fa_FAR)
    global N M epi trail_times

    % Define measurement matrix
    A = get_Psi(N, M, epi);

    % Define sparse vector x
    [Lambda, x] = get_sparse_vector(N * M, 1);
    Lambda_C = setdiff(1:N*M, Lambda);

    % Try "trail_times" times
    results_0 = zeros(trail_times, N * M);
    results_1 = zeros(trail_times, N * M);
    thresholds_0 = [];
    thresholds_1 = [];
    
    noise_sigma = 0.01;
    noise_sigma_2 = noise_sigma ^ 2;
    
    if exist("520BP.mat")
        load("50.mat", "results_0", "results_1");
    else
        % 0 假设
        for T = 1: trail_times
            % y = Ax + n
            n = randn(N, 1) * noise_sigma;
            y_noise = n;
            [x_hat, threshold] = debiased_LASSO(A, y_noise, P_fa_FAR, noise_sigma_2);
            results_0(T, :) = x_hat;
            thresholds_0 = [thresholds_0 threshold];
        end
        % 1 假设
        for T = 1: trail_times
            % y = Ax + n
            n = randn(N, 1) * noise_sigma;
            y_noise = A * x + n;
            [x_hat, threshold] = debiased_LASSO(A, y_noise, P_fa_FAR, noise_sigma_2);
            results_1(T, :) = x_hat;
            thresholds_1 = [thresholds_1 threshold];
        end
    end

    % Distribute of H_0 and H_1
    H_00_distribute = zeros(1, length(Lambda_C));
    H_01_distribute = zeros(1, length(Lambda));
    H_10_distribute = zeros(1, length(Lambda_C));
    H_11_distribute = zeros(1, length(Lambda));

    i1 = 1;
    i2 = 1;
    T_00 = [];
    T_01 = [];
    T_10 = [];
    T_11 = [];
    
    for threshold = 1: trail_times
        x_0_hat = results_0(threshold, :);
        x_1_hat = results_1(threshold, :);
        t_00 = 0;
        t_01 = 0;
        t_10 = 0;
        t_11 = 0;

        for i = 1: length(x)
            if ismember(i, Lambda_C)
                t_00 = t_00 + abs(x_0_hat(i));
                t_10 = t_10 + abs(x_1_hat(i));
                H_00_distribute(i1) = x_0_hat(i);
                H_10_distribute(i1) = x_1_hat(i);
                i1 = i1 + 1;
            elseif ismember(i, Lambda)
                t_01 = t_01 + abs(x_0_hat(i));
                t_11 = t_11 + abs(x_1_hat(i));
                H_01_distribute(i2) = x_0_hat(i);
                H_11_distribute(i2) = x_1_hat(i);
                i2 = i2 + 1;
            end
        end
        T_00 = [T_00 t_00];
        T_01 = [T_01 t_01];
        T_10 = [T_10 t_10];
        T_11 = [T_11 t_11];
    end


    % Draw
    figure(1)
    subplot(2, 2, 1);
    title("Freq histogram of H_0");
    histfit(real(H_00_distribute));
    
    subplot(2, 2, 2);
    title("Freq histogram of H_1");
    histfit(real(H_01_distribute));

    subplot(2, 2, 3);
    title("Freq histogram of H_0");
    histfit(real(H_10_distribute));
    
    subplot(2, 2, 4);
    title("Freq histogram of H_1");
    histfit(real(H_11_distribute));
    
    H_00_mean = mean(H_00_distribute);
    H_00_std = std(H_00_distribute);
    H_01_mean = mean(H_01_distribute);
    H_01_std = std(H_01_distribute);
    H_10_mean = mean(H_10_distribute);
    H_10_std = std(H_10_distribute);
    H_11_mean = mean(H_11_distribute);
    H_11_std = std(H_11_distribute);

    fprintf("H_00: mu = %f, std = %f\n", H_00_mean, H_00_std);
    fprintf("H_01: mu = %f, std = %f\n", H_01_mean, H_01_std);
    fprintf("H_10: mu = %f, std = %f\n", H_10_mean, H_10_std);
    fprintf("H_11: mu = %f, std = %f\n", H_11_mean, H_11_std);

    figure(2);
    subplot(2, 2, 1);
    histfit(real(T_00));
    
    subplot(2, 2, 2);
    histfit(real(T_01));

    subplot(2, 2, 3);
    histfit(real(T_10));
    
    subplot(2, 2, 4);
    histfit(real(T_11));

    % P_fa = normcdf(threshold, H_0_mean, H_0_std);
end