function P_fa = PD_simu(threshold)
    global N method trail_times

    % Define measurement matrix
    B = dftmtx(N);
    
    % Define sparse vector z
    [Lambda, z] = get_sparse_vector(N, 1);
    Lambda_C = setdiff(1:N, Lambda);
    results = zeros(trail_times, N);

    % Calcuate test statistics
    P2 = [];

    if exist("PD_recovery_0_debiasedLASSO_500.mat")
        load("PD_recovery_1_debiasedLASSO_500.mat", "results");
    else
        for T = 1: trail_times
            % y = Bz + n
            n = randn(N, 1) * 0.1;
            y_noise = n;
            
            % z_hat = CS(B, y)
            z_hat = recovery(B, y_noise);
    
            % Calculate test statistic
            results(T, :) = z_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
        z_hat = results(t, :);
        T = 0;
        for i = 1: length(z)
            if ismember(i, Lambda_C)
                H_0_distribute(i1) = z_hat(i);
                i1 = i1 + 1;
            elseif ismember(i, Lambda)
                T = T + abs(z_hat(i));
                H_1_distribute(i2) = z_hat(i);
                i2 = i2 + 1;
            end
        end
        Ts = [Ts T];
    end

    % Draw
    figure(2)
    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