function P_fa = FAR_simu(threshold)
    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
    NN = N;
    results = zeros(trail_times, N * M);

    if exist("FAR_recovery_1_LASSO_5.mat")
        load("FAR_recovery_0_debiasedLASSO_500.mat", "results");
    else
        for T = 1: trail_times
            % y = Ax + n
            n = randn(NN, 1) * 0.1;
            y_noise = A * x + n;
            
            % x_hat = CS(A, y)
            x_hat = recovery(A, y_noise);
            results(T, :) = x_hat;
        end
    end

    % Distribute of H_0 and H_1
    H_0_distribute = zeros(1, length(Lambda_C));
    H_1_distribute = zeros(1, length(Lambda));

    i1 = 1;
    i2 = 1;
    Ts = [];
    
    for t = 1: trail_times
        x_hat = results(t, :);
        T = 0;
        for i = 1: length(x)
            if ismember(i, Lambda_C)
                H_0_distribute(i1) = x_hat(i);
                i1 = i1 + 1;
            elseif ismember(i, Lambda)
                T = T + abs(x_hat(i));
                H_1_distribute(i2) = x_hat(i);
                i2 = i2 + 1;
            end
        end
        Ts = [Ts T];
    end


    % Draw
    figure(1)
    subplot(2, 1, 1);
    title("Freq histogram of H_0");
    histfit(real(H_0_distribute));
    
    subplot(2, 1, 2);
    title("Freq histogram of H_1");
    histfit(real(H_1_distribute));
    
    H_0_mean = mean(H_0_distribute);
    H_0_std = std(H_0_distribute);
    H_1_mean = mean(H_1_distribute);
    H_1_std = std(H_1_distribute);

    fprintf("H_0: mu = %f, std = %f\n", H_0_mean, H_0_std);
    fprintf("H_1: mu = %f, std = %f\n", H_1_mean, H_1_std);
    % P_fa = normcdf(threshold, H_0_mean, H_0_std);
end