Perfomance Comaprison

Perfomance Comaprison#

We consider antipodal signals where

  • \( s_0 = -b \)

  • \( s_1 = b \)

with the a priori probabilities

  • \( \pi_0 = 0.8 \)

  • \( \pi_1 = 1 - \pi_0 = 0.2 \)

We assume \( b > 0 \).

The noise \( \mathbf{n} \) is modeled as a zero-mean Gaussian random variable with variance \( \sigma^2 \), with the probability density function (pdf):

\[ p_{\mathbf{n}}(n) = \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left( -\frac{n^2}{2\sigma^2} \right) \]

In this example, we set

  • \( b = 1 \)

  • \( \sigma = 0.5 \)

import numpy as np
from scipy.special import erf
from scipy.optimize import fsolve
from scipy.integrate import quad
from scipy.special import erfcinv
import pandas as pd

# Given parameters
b = 1
sigma = 0.5 # std, thus, variance sigma^2 = 0.25
pi0 = 0.8
pi1 = 0.2

# Number of simulations
n_samples = 1000000

# Generating s0 and s1
s0 = -b
s1 = b

# Generating noise
noise = np.random.normal(0, np.sqrt(sigma**2), n_samples)

# Simulating received y under each hypothesis
y_given_s0 = s0 + noise
y_given_s1 = s1 + noise

Log-Likelihood Ratio (LLR)#

The likelihood ratio function, which compares the likelihoods of receiving a specific observation \( y \) under two competing hypotheses \( H_1 \) and \( H_0 \), is:

\[ L(y) = \frac{p(y | H_1)}{p(y | H_0)} \]

where:

  • \( p(y | H_1) \) is the probability density function (PDF) of \( y \) given that \( H_1 \) is true.

  • \( p(y | H_0) \) is the PDF of \( y \) given that \( H_0 \) is true.

Given the observation \( y \) in our example, the likelihood functions under the two hypotheses are:

\[ p(y|H_1) = \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left( -\frac{(y - b)^2}{2\sigma^2} \right) \]
\[ p(y|H_0) = \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left( -\frac{(y + b)^2}{2\sigma^2} \right) \]

The Log-Likelihood Ratio (LLR) is simply the logarithm of this likelihood ratio:

\[ \text{LLR}(y) = \log\left(\frac{p(y | H_1)}{p(y | H_0)}\right) \]

The \( \log \) in the LLR typically refers to the natural logarithm, i.e., logarithm to the base \( e \), where \( e \approx 2.71828 \).

So, the expression for the LLR implies:

\[ \text{LLR}(y) = \ln\left(\frac{p(y | H_1)}{p(y | H_0)}\right) \]

where \( \ln \) denotes the natural logarithm.

MAP Criterion#

The Maximum A Posteriori (MAP) criterion is a decision rule used in hypothesis testing that selects the hypothesis with the highest posterior probability.

It incorporates both the likelihood of the observation and the prior probabilities of the hypotheses.

The MAP criterion aims to minimize the overall probability of error by factoring in the prior knowledge about the hypotheses.

Threshold \( \tau_{\text{MAP}} \)#

The threshold \( \tau_{\text{MAP}} \) in the MAP criterion is defined as the ratio of the prior probabilities of the hypotheses \( H_0 \) and \( H_1 \):

\[ \tau_{\text{MAP}} = \frac{\pi_0}{\pi_1} \]

In this example, the threshold is:

\[ \tau_{\text{MAP}} = \frac{0.8}{0.2} = 4 \]

Decision Rule Based on \( \tau_{\text{MAP}} \)#

The likelihood ratio function \( L(y) \), which compares the likelihood of receiving the observation \( y \) under both hypotheses, is given by:

\[ L(y) = \frac{p(y | H_1)}{p(y | H_0)} = \exp\left( \frac{2yb}{\sigma^2} \right) \]

The decision rule based on \( \tau_{\text{MAP}} \) is:

\[\begin{split} L(y) \begin{array}{c} H_1 \\ \gtrless \\ H_0 \end{array} \tau_{\text{MAP}} \end{split}\]

in an alternative presenation:

\[\begin{split} \delta_{\text{MAP}}(y) = \begin{cases} 1, & \text{if } L(y) \geq \tau_{\text{MAP}} \\ 0, & \text{if } L(y) < \tau_{\text{MAP}} \end{cases} \end{split}\]

which means:

  • Decide \( H_1 \) if \( L(y) \geq \tau_{\text{MAP}} \)

  • Decide \( H_0 \) if \( L(y) < \tau_{\text{MAP}} \)

Threshold \( \tau'_{\text{MAP}} \)#

Taking the natural logarithm on both sides:

\[ \frac{2yb}{\sigma^2} \geq \ln 4 \iff y \geq \underbrace{\frac{\sigma^2 \ln 4}{2b}}_{\triangleq \tau'_{\text{MAP}}} \]

\( \tau'_{\text{MAP}} \) is the threshold after taking the natural logarithm of the likelihood ratio:

\[ \tau'_{\text{MAP}} = \ln \left( \frac{\pi_0}{\pi_1} \right) \]

In this example, for \( \pi_0 = 0.8 \) and \( \pi_1 = 0.2 \):

\[ \tau'_{\text{MAP}} = \frac{0.25 \ln 4}{2 \cdot 1} = 0.1733 \]

Decision Rule Based on \( \tau'_{\text{MAP}} \)#

The decision rule based on \( \tau'_{\text{MAP}} \) is:

\[\begin{split} y \begin{array}{c} H_1 \\ \gtrless \\ H_0 \end{array} \tau'_{\text{MAP}} \end{split}\]

in an alternative presenation:

\[\begin{split} \delta'_{\text{MAP}}(y) = \begin{cases} 1, & \text{if } y \geq \tau'_{\text{MAP}} \\ 0, & \text{if } y < \tau'_{\text{MAP}} \end{cases} \end{split}\]

which means:

  • Decide \( H_1 \) if \( y \geq \tau'_{\text{MAP}} \)

  • Decide \( H_0 \) if \( y < \tau'_{\text{MAP}} \)

# MAP criterion threshold
tau_MAP = np.log(pi0 / pi1)

# Decision rule based on MAP criterion
threshold_MAP = (sigma**2 / (2 * b)) * tau_MAP

# Decisions
decisions_s0_MAP = y_given_s0 >= threshold_MAP  # Decisions when H0 is true (y_given_s0)
decisions_s1_MAP = y_given_s1 >= threshold_MAP  # Decisions when H1 is true (y_given_s1)

# Estimating transmitted signals (s_hat) based on decisions
s0_hat_MAP = np.where(decisions_s0_MAP, s1, s0)  # Decided s1 if decision is True, else s0 (under H0)
s1_hat_MAP = np.where(decisions_s1_MAP, s1, s0)  # Decided s1 if decision is True, else s0 (under H1)

# Calculating empirical probabilities
P00_MAP = np.mean(s0_hat_MAP == s0)  # Probability of deciding H0 when H0 is true
P01_MAP = np.mean(s1_hat_MAP == s0)  # Probability of deciding H0 when H1 is true
P10_MAP = np.mean(s0_hat_MAP == s1)  # Probability of deciding H1 when H0 is true
P11_MAP = np.mean(s1_hat_MAP == s1)  # Probability of deciding H1 when H1 is true

# Print results
print(f"MAP decision threshold (threshold_MAP): {threshold_MAP:.4f}")
print(f"MAP, P00 (Decide H0 | H0 true): {P00_MAP:.4f}")
print(f"MAP, P01 (Decide H0 | H1 true): {P01_MAP:.4f}")
print(f"MAP, P10 (Decide H1 | H0 true): {P10_MAP:.4f}")
print(f"MAP, P11 (Decide H1 | H1 true): {P11_MAP:.4f}")

MAP decision threshold (threshold_MAP): 0.1733
MAP, P00 (Decide H0 | H0 true): 0.9904
MAP, P01 (Decide H0 | H1 true): 0.0491
MAP, P10 (Decide H1 | H0 true): 0.0096
MAP, P11 (Decide H1 | H1 true): 0.9509

Bayes#

Bayes Criterion#

The Bayes criterion is a decision-making framework in hypothesis testing that accounts for the cost associated with making incorrect decisions.

By incorporating these costs, the criterion aims to minimize the expected risk, called the Bayes risk.

Cost \( C_{ij} \)#

Let \( C_{ij} \) represent the cost of deciding hypothesis \( D_i \) when the true hypothesis is \( H_j \).

These costs reflect the consequences of making correct/incorrect decisions:

  • \( C_{00} \): Cost of deciding \( D_0 \) when \( H_0 \) is true (correct decision).

  • \( C_{10} \): Cost of deciding \( D_1 \) when \( H_0 \) is true (incorrect decision).

  • \( C_{01} \): Cost of deciding \( D_0 \) when \( H_1 \) is true (incorrect decision).

  • \( C_{11} \): Cost of deciding \( D_1 \) when \( H_1 \) is true (correct decision).

Bayes Risk, \(r\)#

The Bayes risk \( r \) is the expected cost of a decision, accounting for both the prior probabilities of the hypotheses and the decision costs.

It is expressed as:

\[ r = [P_{00}C_{00} + P_{10}C_{10}]\pi_0 + [P_{01}C_{01} + P_{11}C_{11}]\pi_1 \]

where \( P_{ij} \) represents the probability of deciding \( D_i \) (accepting \( H_i \)) when \( H_j \) is true.

Our goal is to minimize the Bayes risk \( r \) by optimizing the decision rule.

Bayes Decision Threshold#

The Bayes decision threshold \( \tau_B \) is derived by balancing the costs and the prior probabilities of the hypotheses.

It is given by:

\[ \tau_B \triangleq \frac{\pi_0 (C_{10} - C_{00})}{\pi_1 (C_{01} - C_{11})} \]

In this example, the costs are assigned as:

  • \( C_{00} = C_{11} = 0 \) (no cost for correct decisions).

  • \( C_{10} = 1 \) (cost of deciding \( D_1 \) when \( H_0 \) is true).

  • \( C_{01} = 2 \) (cost of deciding \( D_0 \) when \( H_1 \) is true).

\( \tau_B \) is calculated as:

\[ \tau_B = \frac{0.8 (1 - 0)}{0.2 (2 - 0)} = 2 \]

Bayes Decision Rule#

The Bayes decision rule, based on the observation \( y \), we compare the likelihood ratio \( L(y) = \frac{p(y|H_1)}{p(y|H_0)} \) to \( \tau_B \):

\[ L(y) \underset{D_0}{\overset{D_1}{\gtrless}} \tau_B \]

alternative presentation

\[\begin{split} \delta_B(y) = \begin{cases} 1, & \text{if } L(y) \geq \tau_B \\ 0, & \text{if } L(y) < \tau_B \end{cases} \end{split}\]

which means:

  • Decide \( H_1 \) if \( L(y) \geq \tau_B \)

  • Decide \( H_0 \) if \( L(y) < \tau_B \)

Threshold \( \tau'_B \)#

Using the log-likelihood ratio, the logarithmic threshold \( \tau'_B \) is:

\[ \tau' = \ln 2 \frac{\sigma^2}{2b} = 0.0866 \]

The decision rule based on \( \tau' \) is:

\[ y \underset{D_0}{\overset{D_1}{\gtrless}} \tau'_B \]

alternatively,

\[\begin{split} \delta'_B (y) = \begin{cases} 1, & \text{if } y \geq \tau'_B \\ 0, & \text{if } y < \tau'_B \end{cases} \end{split}\]

which means:

  • Decide \( H_1 \) if \( y \geq \tau'_B \)

  • Decide \( H_0 \) if \( y < \tau'_B \)

The decision regions are:

  • \( R_1 = \{ y : y \geq \ln 2 \frac{\sigma^2}{2b} \} \) (decide \( H_1 \))

  • \( R_0 = \{ y : y < \ln 2 \frac{\sigma^2}{2b} \} \) (decide \( H_0 \))

# Costs assigned to decisions
C00 = 0  # Cost of deciding D0 when H0 is true (correct decision, no cost)
C11 = 0  # Cost of deciding D1 when H1 is true (correct decision, no cost)
C10 = 1  # Cost of deciding D1 when H0 is true (false positive)
C01 = 2  # Cost of deciding D0 when H1 is true (false negative)

# Bayes decision threshold tau_B calculation
tau_B = (pi0 * (C10 - C00)) / (pi1 * (C01 - C11))

# Bayes decision threshold tau_prime using log-likelihood ratio
tau_prime_B = np.log(tau_B) * (sigma**2 / (2 * b))

# Decisions based on the Bayes threshold
decisions_s0_B = y_given_s0 >= tau_prime_B  # Decide s1 if y >= threshold
decisions_s1_B = y_given_s1 >= tau_prime_B  # Decide s1 if y >= threshold

# Estimating transmitted signals (s_hat) based on decisions for Bayes criterion
s0_hat_B = np.where(decisions_s0_B, s1, s0)  # Decided s1 if decision is True, else s0 (under H0)
s1_hat_B = np.where(decisions_s1_B, s1, s0)  # Decided s1 if decision is True, else s0 (under H1)

# Calculating empirical probabilities for Bayes criterion
P00_Bayes = np.mean(s0_hat_B == s0)  # Probability of deciding H0 when H0 is true
P01_Bayes = np.mean(s1_hat_B == s0)  # Probability of deciding H0 when H1 is true
P10_Bayes = np.mean(s0_hat_B == s1)  # Probability of deciding H1 when H0 is true
P11_Bayes = np.mean(s1_hat_B == s1)  # Probability of deciding H1 when H1 is true

# Print results
print(f"Bayes decision threshold (tau_prime_B): {tau_prime_B:.4f}")
print(f"Bayes, P00 (Decide H0 | H0 true) using Bayes decision threshold: {P00_Bayes:.4f}")
print(f"Bayes, P01 (Decide H0 | H1 true) using Bayes decision threshold: {P01_Bayes:.4f}")
print(f"Bayes, P10 (Decide H1 | H0 true) using Bayes decision threshold: {P10_Bayes:.4f}")
print(f"Bayes, P11 (Decide H1 | H1 true) using Bayes decision threshold: {P11_Bayes:.4f}")

Bayes decision threshold (tau_prime_B): 0.0866
Bayes, P00 (Decide H0 | H0 true) using Bayes decision threshold: 0.9849
Bayes, P01 (Decide H0 | H1 true) using Bayes decision threshold: 0.0338
Bayes, P10 (Decide H1 | H0 true) using Bayes decision threshold: 0.0151
Bayes, P11 (Decide H1 | H1 true) using Bayes decision threshold: 0.9662

Minimax Criterion#

The Minimax criterion is a decision-making approach used in hypothesis testing when the prior probabilities of the hypotheses are unknown or uncertain.

The goal of the minimax strategy is to minimize the maximum possible risk, ensuring the best worst-case performance.

Actual Risk#

When the prior probability is unknown, we must guess it, e.g., we guess the prior probability \( \pi_0 \) is a specific value \( a, 0 \leq a \leq 1 \).

For such a given guess of \( \pi_0 = a \), the actual risk \( r(\pi_0, a) \) is the expected risk based on this guessed prior, is defined as:

\[ r(\pi_0, a) = \pi_0 r_0(a) + (1 - \pi_0) r_1(a) \]

where:

  • \( r_0(a) \) and \( r_1(a) \) are the conditional risks for hypotheses \( H_0 \) and \( H_1 \), respectively, when the decision threshold is controlled by the guessed prior \( a \).

The conditional risk for hypothesis \( H_j \) is given by:

\[ r_j(a) = C_{1j}P_{1j} + C_{0j}P_{0j}, \quad j = 0, 1 \]

Discussion. If the guessed prior \( a \) matches the true prior \( \pi_0 \), then the actual risk is identical to the Bayes risk. That is why the Bayes risk curve intersects with the actual risk curve.

Endpoint Minimax Solution#

In certain cases, the minimax solution corresponds to always choosing a single hypothesis, regardless of the observation.

This situation occurs when the risk associated with one endpoint of the prior probability space is higher than the other.

  • If the minimax risk corresponds to \( \pi_0 = 0 \): The minimax decision is to always choose \( D_1 \) (the hypothesis \( H_1 \)), since the risk associated with \( C_{11} \) is considered the lowest.

    • Decision rule:

      \[ \delta_{MM} (y) = D_1, \quad\forall y \]

    which means: Choose \( D_1 \) for all observations \( y \).

  • If the minimax risk corresponds to \( \pi_0 = 1 \): The minimax decision is to always choose \( D_0 \) (the hypothesis \( H_0 \)), since the risk associated with \( C_{00} \) is minimized.

    • Decision:

    \[ \delta_{MM} (y) = D_0, \quad\forall y \]

    Choose \( D_0 \) for all observations \( y \).

In this approach, we compare the risks at the endpoints (\( \pi_0 = 0 \) and \( \pi_0 = 1 \)) and pick the maximum of \( C_{00} \) or \( C_{11} \) as the minimax risk if the endpoint minimax solution applies.

Decision Threhold based on Differentiable Interior Minimax Solution#

In cases where the minimax solution lies within the interior of the prior probability space (i.e., \( 0 < \pi_0 < 1 \)), we solve for the minimax solution by finding the point where the slope of the actual risk is zero.

Find the optimal ``guessed’’ value \( a_{opt} \) theoretically (no longer a guessing). This corresponds to the condition:

\[ r_0(a_{opt}) = r_1(a_{opt}) \]

In other words, the minimax solution occurs when the conditional risks for both hypotheses are equal.

Find minimax threshold \( \tau_{MM} \) using the obtained \( a_{opt} \):

\[ \tau_{MM} = \frac{a_{opt}(C_{10} - C_{00})}{(1 - a_{opt})(C_{01} - C_{11})} \]

Decision Rule based on \( \tau_{MM} \)#

The MiniMax decision rule, based on the observation \( y \), we compare the likelihood ratio \( L(y) = \frac{p(y|H_1)}{p(y|H_0)} \) to \( \tau_{MM} \):

\[ L(y) \underset{D_0}{\overset{D_1}{\gtrless}} \tau_{MM} \]

alternative presentation

\[\begin{split} \delta_B(y) = \begin{cases} 1, & \text{if } L(y) \geq \tau_{MM} \\ 0, & \text{if } L(y) < \tau_{MM} \end{cases} \end{split}\]

which means:

  • Decide \( H_1 \) if \( L(y) \geq \tau_{MM} \)

  • Decide \( H_0 \) if \( L(y) < \tau_{MM} \)

Decision Rule based on \( \tau'_{MM} \)#

\( \tau'_{MM} \) is the root of this

\[ r_0(\tau'_{MM}) = r_1(\tau'_{MM}) \]

Using the observation \( y \) directly (not the LR \( L(y) \), the decision rule

\[ y \underset{D_0}{\overset{D_1}{\gtrless}} \tau'_{MM} \]

alternative presentation

\[\begin{split} \delta_B(y) = \begin{cases} 1, & \text{if } y \geq \tau'_{MM} \\ 0, & \text{if } y < \tau'_{MM} \end{cases} \end{split}\]

which means:

  • Decide \( H_1 \) if \( y \geq \tau'_{MM} \)

  • Decide \( H_0 \) if \( y < \tau'_{MM} \)

In this exmaple, to find \(\tau'_{MM}\), we solve

\[ \begin{align}\begin{aligned} \pi_0 (\tau'_{MM}) = \pi_1 (\tau'_{MM})\\\iff C_{10}P_{10} + C_{00}P_{00} = C_{11}P_{11} + C_{01}P_{01}\\\iff \text{erf} \left( \frac{b - \tau_{MM}'}{\sqrt{2\sigma}} \right) - \frac{1}{2} \text{erf} \left( \frac{\tau_{MM}' + b}{\sqrt{2\sigma}} \right) = \frac{1}{2} \end{aligned}\end{align} \]

and we obtain \( \tau_{MM}' \approx -0.073 \)

# Define the function f for the Minimax criterion
def f(t_MM_prime):
    term1 = erf((b - t_MM_prime) / np.sqrt(2 * sigma**2))
    term2 = 0.5 * erf((t_MM_prime + b) / np.sqrt(2 * sigma**2))
    return term1 - term2 - 0.5

# Solve for the minimax threshold t_MM_prime using fsolve
t_MM_prime = fsolve(f, 0)[0]  # Initial guess is 0

# Decisions based on the Minimax threshold
decisions_s0_MM = y_given_s0 >= t_MM_prime  # Decide s1 if y >= t_MM_prime (for H0)
decisions_s1_MM = y_given_s1 >= t_MM_prime  # Decide s1 if y >= t_MM_prime (for H1)

# Estimating transmitted signals (s_hat) based on decisions for Minimax criterion
s0_hat_MM = np.where(decisions_s0_MM, s1, s0)  # Decided s1 if decision is True, else s0 (under H0)
s1_hat_MM = np.where(decisions_s1_MM, s1, s0)  # Decided s1 if decision is True, else s0 (under H1)

# Calculating empirical probabilities for Minimax criterion
P00_MM = np.mean(s0_hat_MM == s0)  # Probability of deciding H0 when H0 is true
P01_MM = np.mean(s1_hat_MM == s0)  # Probability of deciding H0 when H1 is true
P10_MM = np.mean(s0_hat_MM == s1)  # Probability of deciding H1 when H0 is true
P11_MM = np.mean(s1_hat_MM == s1)  # Probability of deciding H1 when H1 is true

# Print the results
print(f"Minimax Criterion Threshold (t_MM_prime): {t_MM_prime:.4f}")
print(f"Minimax Criterion, P00 (Decide H0 | H0 true): {P00_MM:.4f}")
print(f"Minimax Criterion, P01 (Decide H0 | H1 true): {P01_MM:.4f}")
print(f"Minimax Criterion, P10 (Decide H1 | H0 true): {P10_MM:.4f}")
print(f"Minimax Criterion, P11 (Decide H1 | H1 true): {P11_MM:.4f}")

Minimax Criterion Threshold (t_MM_prime): -0.0730
Minimax Criterion, P00 (Decide H0 | H0 true): 0.9679
Minimax Criterion, P01 (Decide H0 | H1 true): 0.0158
Minimax Criterion, P10 (Decide H1 | H0 true): 0.0321
Minimax Criterion, P11 (Decide H1 | H1 true): 0.9842

Neyman-Pearson (NP) Criterion#

The Neyman-Pearson (NP) criterion is a framework used in hypothesis testing when one seeks to maximize the probability of detection, \( P_d \), subject to a constraint on the probability of false alarm, \( P_f \).

This is often used in signal detection problems to ensure a desired level of performance in detecting a signal while limiting false detections.

NP Threshold \(\tau_{NP}\)#

We aim to find a threshold \( \tau_{NP} \) that maximizes \( P_d = \Pr(D_1|H_1) \) while satisfying the constraint \( P_f = \Pr(D_1|H_0) \leq \alpha_f \).

To approach this problem, we use a Lagrange multiplier \( \mu \geq 0\) to formulate a function \( F \) that incorporates both the objective of maximizing \( P_d \) and the constraint on \( P_f \):

\[ F = P_d - \mu [P_f - \alpha'] = \text{func} (\mu, \alpha') \]

where \(P_f = \alpha' \leq \alpha_f\).

This inequality leads to a likelihood-ratio test (LRT), where we compare the likelihood ratio \( L(y) \) to a threshold \( \mu \):

\[ L(y) = \frac{p(y|H_1)}{p(y|H_0)} \geq \mu \]

Optimal threshold \(\tau_{NP}\)#

The optimal threshold \( \mu^* \) is chosen such that:

\[ \tau_{NP} = \mu^* \quad \text{where} \quad P_f = \alpha' =: \alpha_f \]

This threshold \( \tau_{NP} \) maximizes \( P_d \) while satisfying the false-alarm constraint.

Neyman-Pearson Decision Rule based on \(\tau_{NP}\)#

The decision rule for the NP criterion can be summarized as:

\[\begin{split} \delta_{NP}(y) = \begin{cases} 1, & L(y) > \tau_{NP} \\ \eta, & L(y) = \tau_{NP} \\ 0, & L(y) < \tau_{NP} \end{cases} \end{split}\]

  • If the likelihood ratio \( L(y) > \tau_{NP} \), we decide \( D_1 \) (i.e., choose hypothesis \( H_1 \)).

  • If \( L(y) < \tau_{NP} \), we decide \( D_0 \) (i.e., choose hypothesis \( H_0 \)).

  • If \( L(y) = \tau_{NP} \), we make a probabilistic decision using a randomization parameter \( \eta \), where the probability of deciding \( D_1 \) is \( \eta \), and the probability of deciding \( D_0 \) is \( 1 - \eta \).

NP Threshold \(\tau'_{NP}\)#

Instead of iterating over \( \alpha'_f \) or solving the Lagrangian function to find the optimal \( \mu \), we directly set the threshold \( \tau_{NP} \) such that \( P_f = \alpha_f \).

Assuming \(P_f = \text{func} (\tau'_{NP})\), we have \(\tau'_{NP} = \text{func}^{-1} (\alpha_f)\).

In this example,

\[ \tau'_{NP} = -b + \sqrt{2\sigma^2} \, \text{erfc}^{-1}(2\alpha_f) \approx -0.36007 \]

Neyman-Pearson Decision Rule based on \(\tau'_{NP}\)#

The NP decision rule is

\[\begin{split} \delta_{NP}(y) = \begin{cases} 1, & y \geq \tau'_{NP} \\ 0, & y < \tau'_{NP} \end{cases} \end{split}\]

which means:

  • Decide \( H_1 \) if \( y \geq \tau'_{NP} \)

  • Decide \( H_0 \) if \( y < \tau_{NP} \)

alpha_f = 0.1  # Desired false alarm probability

# Mathematical expressions for P_f and P_d
def P_f_math(tau_NP):
    func = lambda y: np.exp(-(y + b)**2 / (2 * sigma**2)) / np.sqrt(2 * np.pi * sigma**2)
    return quad(func, tau_NP, np.inf)[0]

def P_d_math(tau_NP):
    func = lambda y: np.exp(-(y - b)**2 / (2 * sigma**2)) / np.sqrt(2 * np.pi * sigma**2)
    return quad(func, tau_NP, np.inf)[0]

# Calculate tau'_NP threshold for Neyman-Pearson criterion
def tau_NP_theoretical(alpha_f, b, sigma):
    return -b + np.sqrt(2 * sigma**2) * erfcinv(2 * alpha_f)

# Calculate tau'_NP
tau_NP = tau_NP_theoretical(alpha_f, b, sigma)

# Decisions based on the Neyman-Pearson threshold
decisions_s0_NP = y_given_s0 >= tau_NP  # Decide H1 if y >= tau'_NP (for H0)
decisions_s1_NP = y_given_s1 >= tau_NP  # Decide H1 if y >= tau'_NP (for H1)

# Estimating transmitted signals (s_hat) based on decisions for Neyman-Pearson criterion
s0_hat_NP = np.where(decisions_s0_NP, s1, s0)  # Decided s1 if decision is True, else s0 (under H0)
s1_hat_NP = np.where(decisions_s1_NP, s1, s0)  # Decided s1 if decision is True, else s0 (under H1)

# Calculating empirical probabilities for Neyman-Pearson criterion
P00_NP = np.mean(s0_hat_NP == s0)  # Probability of deciding H0 when H0 is true
P01_NP = np.mean(s1_hat_NP == s0)  # Probability of deciding H0 when H1 is true
P10_NP = np.mean(s0_hat_NP == s1)  # Probability of deciding H1 when H0 is true (false alarm)
P11_NP = np.mean(s1_hat_NP == s1)  # Probability of deciding H1 when H1 is true (detection)

# Print the results
print(f"Neyman-Pearson threshold tau'_NP = {tau_NP}")
print(f"Neyman-Pearson Criterion, P00 (Decide H0 | H0 true): {P00_NP:.4f}")
print(f"Neyman-Pearson Criterion, P01 (Decide H0 | H1 true): {P01_NP:.4f}")
print(f"Neyman-Pearson Criterion, P10 (Decide H1 | H0 true): {P10_NP:.4f}")
print(f"Neyman-Pearson Criterion, P11 (Decide H1 | H1 true): {P11_NP:.4f}")

Neyman-Pearson threshold tau'_NP = -0.3592242172276997
Neyman-Pearson Criterion, P00 (Decide H0 | H0 true): 0.8998
Neyman-Pearson Criterion, P01 (Decide H0 | H1 true): 0.0033
Neyman-Pearson Criterion, P10 (Decide H1 | H0 true): 0.1002
Neyman-Pearson Criterion, P11 (Decide H1 | H1 true): 0.9967

Error Probability \(P_E\)#

The probability of error \( P_E \) for each criterion (MAP, Bayes, Minimax, and Neyman-Pearson) is defined as:

\[ P_E = \pi_0 P_{10} + \pi_1 P_{01} \]

where:

  • \( \pi_0 \) is the prior probability of \( H_0 \) (usually the probability that \( H_0 \) is true).

  • \( \pi_1 \) is the prior probability of \( H_1 \).

  • \( P_{10} \) is the probability of deciding \( H_1 \) when \( H_0 \) is true (false alarm, \( P_f \)).

  • \( P_{01} \) is the probability of deciding \( H_0 \) when \( H_1 \) is true (missed detection, \( P_m \)).

# Prior probabilities
pi_0 = 0.8  # Probability of H0
pi_1 = 0.2  # Probability of H1

# Probability of error formula
def compute_P_E(pi_0, pi_1, P10, P01):
    return pi_0 * P10 + pi_1 * P01


# For MAP Criterion
P_E_MAP = compute_P_E(pi_0, pi_1, P10_MAP, P01_MAP)

# For Bayes Criterion
P_E_Bayes = compute_P_E(pi_0, pi_1, P10_Bayes, P01_Bayes)

# For Minimax Criterion 
P_E_MM = compute_P_E(pi_0, pi_1, P10_MM, P01_MM)

# For Neyman-Pearson Criterion
P_E_NP = compute_P_E(pi_0, pi_1, P10_NP, P01_NP)

# Print the results
print(f"Probability of Error for MAP Criterion: {P_E_MAP:.4f}")
print(f"Probability of Error for Bayes Criterion: {P_E_Bayes:.4f}")
print(f"Probability of Error for Minimax Criterion: {P_E_MM:.4f}")
print(f"Probability of Error for Neyman-Pearson Criterion: {P_E_NP:.4f}")

Probability of Error for MAP Criterion: 0.0175
Probability of Error for Bayes Criterion: 0.0188
Probability of Error for Minimax Criterion: 0.0289
Probability of Error for Neyman-Pearson Criterion: 0.0808

Comparison Table#

# Create a dictionary with the metrics
metrics_data = {
    'Criterion': ['MAP', 'Bayes', 'Minimax', 'Neyman-Pearson'],
    'P00': [P00_MAP, P00_Bayes, P00_MM, P00_NP],
    'P01': [P01_MAP, P01_Bayes, P01_MM, P01_NP],
    'P10 (P_f)': [P10_MAP, P10_Bayes, P10_MM, P10_NP],
    'P11 (P_d)': [P11_MAP, P11_Bayes, P11_MM, P11_NP],
    'PE': [P_E_MAP, P_E_Bayes, P_E_MM, P_E_NP]
}

# Create a DataFrame
metrics_df = pd.DataFrame(metrics_data)

# Display the DataFrame
print(metrics_df.to_string(index=False))

     Criterion      P00      P01  P10 (P_f)  P11 (P_d)       PE
           MAP 0.990379 0.049053   0.009621   0.950947 0.017507
         Bayes 0.984940 0.033830   0.015060   0.966170 0.018814
       Minimax 0.967867 0.015800   0.032133   0.984200 0.028866
Neyman-Pearson 0.899811 0.003265   0.100189   0.996735 0.080804