How can I calculate the probability of false detection?

Hello everyone!

I need to justify for my dissertation the problems of false detection of a signal by a normal distribution (Gaussian) and build a graph where I should get a decreasing exponent in the interval for P from 10^-8 to 10^-1, and for alpha squared (alpha^2) from 10 to 100. The sigma dispersion = from 10 to 100. The threshold value of the signal n2 = 1.5, from which the normal value is integrated to infinity to calculate the probability of P.

I wrote the following code:

% Given parameters
sigma2_values = linspace (0.1, 0.01, 100); % dispersion values from 0.01 to 1000
n2 = 1.5; % threshold value
P_loznoe = zeros(length(sigma2_values), 1); % Initialize array for P_false

% Calculate P_false for each value of sigma^2 for fixed n2
for j = 1: length(sigma2_values)
    sigma2 = sigma2_values(j); % use current value of sigma^2
    alpha2 = 1/sigma2; % Calculate alpha^2
    
    % Calculate integral of P(x) from n2 to infinity
    integrand = @(x) (1 / (sqrt(2 * pi * sigma2))) .* exp(-((x.^2) / (2 * sigma2)));
    P_loznoe(j) = integral(integrand, n2, Inf); % Calculate the integral
end

% Calculate alpha^2 for each sigma^2
alpha2_values = 1 ./sigma2_values; % alpha^2 = 1/sigma^2

% Plot P_false vs. alpha^2
figure;
semilogy(alpha2_values, P_loznoe, 'r', 'LineWidth', 2); % Logarithmic scale on the Y axis
title('False discovery rate vs. \alpha^2');
xlabel('\alpha^2');
ylabel('P_{false}');
xlim([10 100]); % Set limits on the X axis
ylim([10^(-8) 10^(-1)]); % Set limits on the Y axis to expand the grid
grid on; % Grid on

But for some reason my graph is not in the specified interval and not in the form of a decreasing exponent, but in the form of a linear decrease.

How can this be fixed?

Thanks in advance!

[Torsten]

You are aware that you chose a logarithmic scale for the y-axis ?

Answer 1

[Sam Chak]

@Suleyman Aliyev, Do you expect the graph to behave like this?

% Given parameters
sigma2_values   = linspace(0.1, 0.01, 100);         % dispersion values ​​from 0.01 to 1000
n2              = 1.5;                              % threshold value
P_loznoe        = zeros(length(sigma2_values), 1);  % Initialize array for P_false

% Calculate P_false for each value of sigma^2 for fixed n2
for j = 1:length(sigma2_values)
    sigma2      = sigma2_values(j);                 % use current value of sigma^2
    alpha2      = 1/sigma2;                         % Calculate alpha^2
    
    % Calculate integral of P(x) from n2 to infinity
    integrand   = @(x) (1/(sqrt(2*pi*sigma2))) .* exp(-((x.^2)/(2*sigma2)));
    P_loznoe(j) = integral(integrand, n2, Inf);     % Calculate the integral
end

% Calculate alpha^2 for each sigma^2
alpha2_values   = 1./sigma2_values;                 % alpha^2 = 1/sigma^2

% Plot P_false vs. alpha^2
figure;
% semilogy(alpha2_values, P_loznoe, 'r', 'LineWidth', 2); % Logarithmic scale on the Y axis
semilogx(alpha2_values, P_loznoe, 'r', 'LineWidth', 2)
title('False discovery rate vs. \alpha^2');
xlabel('\alpha^2');
ylabel('P_{false}');
xlim([10 100]); % Set limits on the X axis
% ylim([10^(-8) 10^(-1)]); % Set limits on the Y axis to expand the grid
grid on; % Grid on

●

[Suleyman Aliyev]

@Sam Chak Yes, that's right. Can you please explain my mistake to me? I don't have much experience on MATLAB, therefore.

[Sam Chak]

Hi @Suleyman Aliyev

In your original code, the y-axis is set to a logarithmic scale. The command semilogy(x, y) plots the x-coordinates using a standard linear scale on the x-axis and the y-coordinates using a base-10 logarithmic scale on the y-axis.

Since you expect the graph to decay exponentially, I suspect that you intended to use the command semilogx(x, y), which plots the x-coordinates using a base-10 logarithmic scale on the x-axis. In fact, the common plot(x, y) function also exhibits the exponential decay pattern.

[Suleyman Aliyev]

@Sam Chak thanks a lot for your explanations. You've been very helpful.

[Sam Chak]

You're welcome, @Suleyman Aliyev. If you find the explanation helpful, please consider clicking 'Accept' ✔ on the answer and voting 👍 for it. Your support is greatly appreciated!

[Suleyman Aliyev]

@Sam Chak Please tell me if I have a question.:

When you add ylim([10^(-8) 10^(-1)]);, the graph is not being plotted, but I need the graph for P in Y to have values not from 0 to 1.2*10^-6, but from 10^-12 to 10^-1, because what this code builds to 0 is not very suitable, since this is impossible. this is the perfect signal being received.

[Sam Chak]

Hi @Suleyman Aliyev,

ylim([10^(-8) 10^(-1)]); % Set limits on the Y axis to expand the grid

I did not include this line to set the limits on the Y-axis, as it was originally present in your code. However, I disabled it in my response to better display the exponential decay pattern. If you wish to set the limits from to on the Y axis, please feel free to do so.

ylim([1e-12 1e-1])

[Suleyman Aliyev]

@Sam Chak It doesn't work.

[Sam Chak]

@Suleyman Aliyev, Hmm... This line sets the limits on the Y-axis for display purposes, as you previously requested. However, it does not affect the behavior of the graph (red curve) in any significant way.

ylim([1e-12 1e-1])

Your original code remains unchanged. If the graph does not appear as expected, you should verify whether the equations are correct according to the textbooks.

[Suleyman Aliyev]

@Sam Chak I need to switch to a logarithmic scale along the ordinate axis and a linear scale along the abscissa axis, where I get approximately a (more extended) graph, as in the code that you provided, is there any way to do this?

[Torsten]

I need to switch to a logarithmic scale along the ordinate axis and a linear scale along the abscissa axis

Then your original code which used "semilogy" was correct.

[Sam Chak]

@Suleyman Aliyev, Technically, your original code using semilogy() accomplishes exactly that.

semilogy(alpha2_values, P_loznoe, 'r', 'LineWidth', 2)

Additionally, the definite integral can be evaluated analytically.

.

You may be able to simplify the code without performing numerical integration in loops. Also, consider examining the error function, erf(), as it should theoretically provide a more accurate result.

[Suleyman Aliyev]

@Sam Chak For 10, the probability is 1x10 -6. For 20, the probability should drop by at least two orders of magnitude, i.e. 10-8. I need all this to be shown, and not as for 20, again the same order is 0.2x10- 6. Is there any way to do this? Maybe the problem is in the graph display, which is why there is a misunderstanding?

[Torsten]

Your original graph shows a probability of 1e-6 for 10 and appr. 1e-11 for 20. So I don't understand what your problem is.

[Suleyman Aliyev]

@Torsten this is not marked on the graph, I would like it to appear on the graph, so that the graph is more "expanded". Otherwise, according to the schedule, it seems that he is already at zero.

[Suleyman Aliyev]

I made it so that the dots show what is not 0 at 20, but I would like it to be shown more clearly on the graph so that it would be clear even without these dots. Is it possible to do this?

[Torsten]

I don't understand why you write "I need to switch to a logarithmic scale along the ordinate axis and a linear scale along the abscissa axis" and now refer to the graph where the abszissa has logarithmic scale and the ordinate has linear scale.

The plot you have to use is your "semilogy" plot right at the beginning from which you can easily identify the small probability values for larger abzissa values.

[Suleyman Aliyev]

Everything worked out! Thank you all very much for your help!