Probabilistic Models for Physical Layer Security: A Complete Guide for AI Cybersecurity Engineers

What is Physical Layer Security?

Most security engineers spend their careers thinking about encryption keys, certificates, and cryptographic protocols. The assumption is always the same — you have data, you scramble it mathematically, and anyone without the key gets nonsense. That is application layer security, and for decades it has been the dominant paradigm. This article is part of the Scientias AI Labs research hub on Probabilistic Control Engineering for Generative AI.

But there is a fundamentally different way to think about security — one that does not rely on mathematical puzzles or secret keys at all. Physical Layer Security uses the actual physics of wireless communication to make eavesdropping impossible. Not difficult. Not computationally expensive. Mathematically impossible.

The intuition is this. Every wireless channel between two points is physically unique. The channel between you and your intended receiver has specific multipath components, specific noise characteristics, and specific fading patterns that are different from the channel between you and any eavesdropper. If your legitimate receiver consistently experiences a better channel than the eavesdropper — even on average — then information theory guarantees that you can transmit information the eavesdropper cannot recover, regardless of how powerful their hardware is.

This guarantee does not weaken as computers get faster. It does not break when quantum computers arrive. It is rooted in the physics of electromagnetic propagation, not in the difficulty of factoring large numbers.

Why Traditional Cryptography is Not Enough

Classical cryptography is built on computational hardness assumptions. RSA is secure because factoring large numbers is hard today. Elliptic curve cryptography is secure because the discrete logarithm problem is hard today. The word today is doing a lot of work in those sentences.

Shor’s algorithm, running on a sufficiently powerful quantum computer, breaks RSA and ECC efficiently. Post-quantum cryptography is responding to this threat, but it introduces new assumptions and new complexities. Physical Layer Security provides a defense that is orthogonal to all of this — it does not make assumptions about computational hardness at all.

For wireless systems in particular, PLS is becoming increasingly important as the attack surface expands. Every IoT device, every sensor node, every wireless industrial controller is a potential target. Many of these devices are too resource-constrained to run heavy cryptographic protocols. PLS offers a path to security that can be implemented in hardware with minimal overhead.

How PLS Differs from Application Layer Security

Think of it this way. Application layer security is like putting a lock on a safe. Physical layer security is like making the safe invisible to anyone who is not standing in exactly the right position. One approach depends on how strong the lock is. The other depends on the geometry of the space.

---

The Mathematics Behind PLS

Secrecy Capacity — The Core Concept

python

import numpy as np
import matplotlib.pyplot as plt

def secrecy_capacity(snr_main, snr_eve):
    """
    Calculate secrecy capacity of Gaussian wiretap channel
    
    Args:
        snr_main: SNR at legitimate receiver (linear scale)
        snr_eve: SNR at eavesdropper (linear scale)
    
    Returns:
        Secrecy capacity in bits per channel use
    """
    capacity_main = np.log2(1 + snr_main)
    capacity_eve = np.log2(1 + snr_eve)
    cs = capacity_main - capacity_eve
    return max(0, cs)

# Visualize secrecy capacity vs SNR advantage
snr_eve_dB = 5
snr_eve = 10**(snr_eve_dB/10)
snr_main_dB = np.linspace(0, 30, 100)
snr_main = 10**(snr_main_dB/10)

cs_values = [secrecy_capacity(sm, snr_eve) 
             for sm in snr_main]

plt.figure(figsize=(10, 5))
plt.plot(snr_main_dB, cs_values, 
         linewidth=2, color='steelblue')
plt.axvline(x=snr_eve_dB, color='red', 
            linestyle='--', label=f'Eve SNR = {snr_eve_dB} dB')
plt.xlabel('Main Channel SNR (dB)')
plt.ylabel('Secrecy Capacity (bits/channel use)')
plt.title('Secrecy Capacity vs Main Channel SNR')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()

# Demonstration
for main_dB in [3, 5, 10, 20]:
    main = 10**(main_dB/10)
    cs = secrecy_capacity(main, snr_eve)
    print(f"Main={main_dB}dB, Eve={snr_eve_dB}dB → "
          f"Cs={cs:.3f} bits/use")

What the code does: This function takes two numbers as input — your legitimate receiver’s signal strength and the eavesdropper’s signal strength — and returns the secrecy capacity. If you run the demonstration at the bottom, you will see that when the main channel SNR is below the eavesdropper’s SNR of 5 dB, the secrecy capacity is zero. Once the main channel exceeds 5 dB, secrecy capacity becomes positive and grows as the advantage increases. The plot makes this threshold effect visually clear.

The notation $[x]^+$ means take the maximum of $x$x and zero. Secrecy capacity cannot be negative — if the eavesdropper has the better channel, you simply get zero secure capacity, not negative capacity.

---

Why Probabilistic Models?

Stochastic Channel Modeling

python

import numpy as np
import matplotlib.pyplot as plt

def simulate_fading_channels(n_samples=10000, 
                              sigma_main=1.0, 
                              sigma_eve=0.7):
    """
    Simulate Rayleigh fading channel realizations
    for both legitimate and eavesdropper channels
    
    Args:
        n_samples: Number of channel realizations
        sigma_main: Main channel scale parameter
        sigma_eve: Eavesdropper channel scale parameter
    
    Returns:
        Dictionary of channel samples and secrecy metrics
    """
    # Generate complex Gaussian channel coefficients
    h_main = (np.random.normal(0, sigma_main, n_samples) + 
              1j * np.random.normal(0, sigma_main, n_samples))
    
    h_eve = (np.random.normal(0, sigma_eve, n_samples) + 
             1j * np.random.normal(0, sigma_eve, n_samples))
    
    # Channel gains (squared magnitudes)
    gain_main = np.abs(h_main)**2
    gain_eve = np.abs(h_eve)**2
    
    # Assume transmit SNR of 10 dB
    tx_snr = 10
    snr_main = tx_snr * gain_main
    snr_eve = tx_snr * gain_eve
    
    # Instantaneous secrecy capacity
    cs = np.maximum(0, 
         np.log2(1 + snr_main) - np.log2(1 + snr_eve))
    
    return {
        'gain_main': gain_main,
        'gain_eve': gain_eve,
        'cs': cs,
        'mean_cs': np.mean(cs),
        'sop': np.mean(cs == 0),
        'p_nonzero': np.mean(snr_main > snr_eve)
    }

# Run simulation
results = simulate_fading_channels(n_samples=50000)

fig, axes = plt.subplots(1, 3, figsize=(15, 5))

# Channel gain distributions
axes[0].hist(results['gain_main'], bins=50, 
             density=True, alpha=0.7, label='Main')
axes[0].hist(results['gain_eve'], bins=50, 
             density=True, alpha=0.7, label='Eve')
axes[0].set_xlabel('Channel Gain')
axes[0].set_ylabel('Probability Density')
axes[0].set_title('Channel Gain Distributions')
axes[0].legend()

# Secrecy capacity distribution
axes[1].hist(results['cs'], bins=50, density=True,
             color='steelblue', alpha=0.8)
axes[1].set_xlabel('Secrecy Capacity (bits/use)')
axes[1].set_ylabel('Probability Density')
axes[1].set_title('Secrecy Capacity Distribution')

# Cumulative distribution
cs_sorted = np.sort(results['cs'])
cdf = np.arange(1, len(cs_sorted)+1) / len(cs_sorted)
axes[2].plot(cs_sorted, cdf, linewidth=2)
axes[2].set_xlabel('Secrecy Capacity (bits/use)')
axes[2].set_ylabel('CDF')
axes[2].set_title('CDF of Secrecy Capacity')

plt.tight_layout()
plt.show()

print(f"Mean Secrecy Rate: {results['mean_cs']:.4f} bits/use")
print(f"Secrecy Outage Probability: {results['sop']:.4f}")
print(f"P(nonzero secrecy): {results['p_nonzero']:.4f}")

What the code does: This simulation generates 50,000 random channel realizations for both the legitimate receiver and the eavesdropper. Each realization represents one snapshot of the wireless channel — one moment in time when the channel has a specific fading coefficient. The code then computes the secrecy capacity for each snapshot and builds a statistical picture. The three plots show how channel gains are distributed, how secrecy capacity varies across realizations, and the cumulative distribution function that tells you what fraction of the time secrecy capacity exceeds any given value.

What the math means: Wireless channels are random processes. Every time you transmit, the channel coefficient is drawn from a probability distribution shaped by the physical environment. The key insight for PLS is that we care about the joint behavior of two random channels — the legitimate channel and the eavesdropper’s channel. If the legitimate channel is almost always stronger, then secrecy is almost always available. If the channels are comparable, secrecy becomes intermittent. Probabilistic models let us characterize this statistically and design systems that guarantee secrecy with a specified probability.$$f(h) = \frac{2h}{\Omega}e^{-h^2/\Omega}, \quad h \geq 0$$For Rayleigh fading, the channel magnitude follows this distribution where $\Omega = E[h^2]$ is the average channel power. The exponential decay means large channel gains are rare but possible — and it is precisely during those large gain moments that the highest secrecy rates are achievable.

---

Key Probabilistic Models in PLS

Rayleigh Fading Channel Model

python

import numpy as np
from scipy.stats import rayleigh, nakagami
import matplotlib.pyplot as plt

class ChannelModels:
    """
    Collection of fading channel models for PLS analysis
    """
    
    @staticmethod
    def rayleigh_samples(mean_power, n_samples):
        """
        Generate Rayleigh fading channel gain samples
        
        Args:
            mean_power: Average channel power (Omega)
            n_samples: Number of samples
        
        Returns:
            Channel power gain samples
        """
        sigma = np.sqrt(mean_power / 2)
        h_real = np.random.normal(0, sigma, n_samples)
        h_imag = np.random.normal(0, sigma, n_samples)
        return np.abs(h_real + 1j*h_imag)**2
    
    @staticmethod
    def nakagami_samples(mean_power, m_param, n_samples):
        """
        Generate Nakagami-m fading channel gain samples
        
        Args:
            mean_power: Average channel power
            m_param: Nakagami shape parameter (m>=0.5)
                     m=1 gives Rayleigh, m->inf gives AWGN
            n_samples: Number of samples
        
        Returns:
            Channel power gain samples
        """
        # Gamma distribution gives Nakagami power
        shape = m_param
        scale = mean_power / m_param
        return np.random.gamma(shape, scale, n_samples)
    
    @staticmethod
    def rician_samples(mean_power, K_factor, n_samples):
        """
        Generate Rician fading channel gain samples
        
        Args:
            mean_power: Average channel power
            K_factor: Rician K-factor (ratio LOS/scattered)
                      K=0 gives Rayleigh, K->inf gives AWGN
            n_samples: Number of samples
        
        Returns:
            Channel power gain samples
        """
        # LOS component power
        s = np.sqrt(K_factor * mean_power / (K_factor + 1))
        sigma = np.sqrt(mean_power / (2 * (K_factor + 1)))
        
        h_real = np.random.normal(s, sigma, n_samples)
        h_imag = np.random.normal(0, sigma, n_samples)
        return np.abs(h_real + 1j*h_imag)**2


# Compare models
n_samples = 100000
mean_power = 1.0

rayleigh_gains = ChannelModels.rayleigh_samples(
    mean_power, n_samples)
nakagami_gains = ChannelModels.nakagami_samples(
    mean_power, m_param=2.0, n_samples=n_samples)
rician_gains = ChannelModels.rician_samples(
    mean_power, K_factor=3.0, n_samples=n_samples)

# Compare secrecy capacities under different models
tx_snr = 10
eve_gains = ChannelModels.rayleigh_samples(
    0.5, n_samples)  # Weaker eve channel

models = {
    'Rayleigh': rayleigh_gains,
    'Nakagami-m (m=2)': nakagami_gains,
    'Rician (K=3)': rician_gains
}

plt.figure(figsize=(12, 5))

for i, (name, gains) in enumerate(models.items()):
    snr_main = tx_snr * gains
    snr_eve = tx_snr * eve_gains
    cs = np.maximum(0, 
         np.log2(1 + snr_main) - np.log2(1 + snr_eve))
    
    cs_sorted = np.sort(cs)
    cdf = np.arange(1, len(cs_sorted)+1) / len(cs_sorted)
    plt.plot(cs_sorted, cdf, linewidth=2, label=name)
    
    print(f"{name}:")
    print(f"  Mean Cs: {np.mean(cs):.4f} bits/use")
    print(f"  SOP: {np.mean(cs==0):.4f}")

plt.xlabel('Secrecy Capacity (bits/channel use)')
plt.ylabel('CDF')
plt.title('Secrecy Capacity CDF — Channel Model Comparison')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()

What the code does: This class implements three different fading channel models and lets you compare their impact on secrecy performance. Run it and you will get printed statistics and a CDF plot comparing how secrecy capacity behaves under Rayleigh, Nakagami-m, and Rician fading. The Nakagami-m model with m=2 concentrates channel gains more tightly around the mean, reducing the probability of very deep fades. The Rician model with K=3 has a strong line-of-sight component that keeps channel gains more consistently high. Both improve secrecy performance compared to pure Rayleigh fading.

What the math means: Different physical environments produce different fading statistics. Rayleigh fading happens when there is no dominant propagation path — all the signal components arrive from different directions and add randomly. This is common in urban areas with no line of sight. Nakagami-m generalizes Rayleigh by allowing the fading severity to vary through the parameter m. Rician fading occurs when there is a dominant line-of-sight path alongside scattered components — common in suburban or indoor environments. Choosing the right model for your deployment environment is crucial for accurate PLS analysis.

Nakagami-m PDF:$$f(h) = \frac{2m^m h^{2m-1}}{\Gamma(m)\Omega^m}e^{-mh^2/\Omega}$$Rician PDF:$$f(h) = \frac{2h(K+1)}{\Omega}e^{-K-\frac{(K+1)h^2}{\Omega}}I_0\left(2h\sqrt{\frac{K(K+1)}{\Omega}}\right)$$Where $I_0$​ is the modified Bessel function of the first kind and $K$ is the Rician factor representing the ratio of LOS power to scattered power.

---

Secrecy Metrics — How We Measure PLS

python

import numpy as np
import matplotlib.pyplot as plt

class PLSMetrics:
    """
    Complete suite of Physical Layer Security metrics
    """
    
    def __init__(self, snr_main_samples, snr_eve_samples):
        """
        Args:
            snr_main_samples: Array of main channel SNR samples
            snr_eve_samples: Array of eavesdropper SNR samples
        """
        self.snr_m = snr_main_samples
        self.snr_e = snr_eve_samples
        self.cs = np.maximum(0,
            np.log2(1 + self.snr_m) - 
            np.log2(1 + self.snr_e))
    
    def secrecy_outage_probability(self, target_rate=1.0):
        """P(Cs < Rs) — probability of insecure transmission"""
        return np.mean(self.cs < target_rate)
    
    def average_secrecy_rate(self):
        """E[Cs] — expected secrecy capacity"""
        return np.mean(self.cs)
    
    def probability_nonzero_secrecy(self):
        """P(Cs > 0) — probability secure link exists"""
        return np.mean(self.snr_m > self.snr_e)
    
    def secrecy_throughput(self, target_rate=1.0):
        """
        Rs * P(Cs >= Rs) — effective secure throughput
        """
        return target_rate * (1 - self.secrecy_outage_probability(
            target_rate))
    
    def optimal_secrecy_rate(self, rate_range=None):
        """
        Find rate that maximizes secrecy throughput
        """
        if rate_range is None:
            rate_range = np.linspace(0.1, 5.0, 100)
        
        throughputs = [self.secrecy_throughput(r) 
                      for r in rate_range]
        optimal_idx = np.argmax(throughputs)
        
        return rate_range[optimal_idx], throughputs[optimal_idx]
    
    def full_report(self):
        """Print complete PLS performance report"""
        opt_rate, opt_throughput = self.optimal_secrecy_rate()
        
        print("=" * 50)
        print("Physical Layer Security Performance Report")
        print("=" * 50)
        print(f"Average Secrecy Rate:      {self.average_secrecy_rate():.4f} bits/use")
        print(f"P(Nonzero Secrecy):        {self.probability_nonzero_secrecy():.4f}")
        print(f"SOP at Rs=1.0:             {self.secrecy_outage_probability(1.0):.4f}")
        print(f"SOP at Rs=2.0:             {self.secrecy_outage_probability(2.0):.4f}")
        print(f"Optimal Secrecy Rate:      {opt_rate:.4f} bits/use")
        print(f"Max Secrecy Throughput:    {opt_throughput:.4f} bits/use")
        print("=" * 50)


# Example usage
n_samples = 100000
tx_snr = 10

# Generate channels
h_main = np.random.rayleigh(np.sqrt(tx_snr/2), n_samples)
h_eve = np.random.rayleigh(np.sqrt(tx_snr*0.3/2), n_samples)

snr_main = h_main**2
snr_eve = h_eve**2

metrics = PLSMetrics(snr_main, snr_eve)
metrics.full_report()

# Plot secrecy throughput vs rate
rate_range = np.linspace(0.1, 5.0, 100)
throughputs = [metrics.secrecy_throughput(r) 
               for r in rate_range]

plt.figure(figsize=(10, 5))
plt.plot(rate_range, throughputs, linewidth=2)
plt.xlabel('Target Secrecy Rate (bits/channel use)')
plt.ylabel('Secrecy Throughput (bits/channel use)')
plt.title('Secrecy Throughput vs Target Rate')
opt_rate, opt_tp = metrics.optimal_secrecy_rate()
plt.axvline(x=opt_rate, color='red', linestyle='--',
            label=f'Optimal Rate = {opt_rate:.2f}')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()

What the code does: The PLSMetrics class is a complete toolkit for measuring PLS performance. Feed it two arrays of SNR samples — one for the legitimate channel and one for the eavesdropper — and it computes every metric you need. The full_report() method prints a clean summary. The optimal_secrecy_rate() method finds the transmission rate that maximizes the effective secure throughput — a practically important result because transmitting too fast causes frequent outages while transmitting too slowly wastes capacity.

What the metrics mean: Different PLS metrics answer different engineering questions. Secrecy Outage Probability answers — how often does my secure link fail? Average Secrecy Rate answers — what is my average secure throughput? Probability of Nonzero Secrecy answers — how often does a secure channel even exist? Secrecy Throughput combines rate and reliability into a single number that captures the real engineering tradeoff — transmit faster but risk more outages, or transmit slower but stay reliable.

Secrecy Throughput:$$\tau_s = R_s \cdot \Pr[C_s \geq R_s] = R_s \cdot (1 – P_{out}(R_s))$$This metric captures the fundamental tradeoff in PLS system design. Increasing $R_s$ increases the per-use payload but also increases outage probability. The optimal $R_s^*$​ maximizes $\tau_s$​ and depends on the statistical properties of both channels.

---

Stochastic Geometry for PLS

python

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import exp1

class StochasticGeometryPLS:
    """
    PLS analysis using stochastic geometry
    Models eavesdroppers as Poisson Point Process
    """
    
    def __init__(self, lambda_eve, tx_power=1.0, 
                 noise_power=0.1, path_loss_exp=3.5):
        """
        Args:
            lambda_eve: Density of eavesdroppers (per m^2)
            tx_power: Transmitter power
            noise_power: Noise power
            path_loss_exp: Path loss exponent (2-4 typical)
        """
        self.lambda_e = lambda_eve
        self.P = tx_power
        self.N0 = noise_power
        self.alpha = path_loss_exp
    
    def deploy_network(self, area_radius=100, 
                       n_legitimate=1):
        """
        Deploy random network realization
        
        Returns:
            Dictionary with node positions
        """
        # Eavesdroppers as PPP
        area = np.pi * area_radius**2
        n_eves = np.random.poisson(
            self.lambda_e * area)
        
        angles = np.random.uniform(0, 2*np.pi, n_eves)
        radii = area_radius * np.sqrt(
            np.random.uniform(0, 1, n_eves))
        
        eve_x = radii * np.cos(angles)
        eve_y = radii * np.sin(angles)
        
        return {
            'n_eves': n_eves,
            'eve_positions': np.column_stack([eve_x, eve_y])
        }
    
    def compute_sir(self, node_positions, 
                    target_distance=10.0):
        """
        Compute Signal-to-Interference-plus-Noise Ratio
        for eavesdroppers at given positions
        """
        # Signal power at target receiver
        signal_power = (self.P * 
                       target_distance**(-self.alpha))
        
        # Distances from transmitter to eavesdroppers
        distances = np.sqrt(
            node_positions[:, 0]**2 + 
            node_positions[:, 1]**2)
        distances = np.maximum(distances, 0.1)
        
        # Received power at each eavesdropper
        eve_powers = self.P * distances**(-self.alpha)
        
        return signal_power, eve_powers
    
    def monte_carlo_sop(self, target_rate=1.0,
                        target_distance=10.0,
                        n_trials=1000,
                        area_radius=100):
        """
        Monte Carlo estimation of secrecy outage
        probability with random eavesdropper locations
        """
        outage_count = 0
        
        for _ in range(n_trials):
            # Deploy random network
            network = self.deploy_network(area_radius)
            
            if network['n_eves'] == 0:
                continue
            
            # Generate Rayleigh fading
            h_main = np.random.rayleigh(1.0)
            h_eves = np.random.rayleigh(
                1.0, network['n_eves'])
            
            signal_power, eve_path_losses = self.compute_sir(
                network['eve_positions'], target_distance)
            
            # SNRs
            snr_main = (self.P * h_main**2 * 
                       target_distance**(-self.alpha) / 
                       self.N0)
            
            snr_eves = (h_eves**2 * eve_path_losses / 
                       self.N0)
            
            # Most dangerous eavesdropper
            snr_best_eve = np.max(snr_eves)
            
            # Secrecy capacity
            cs = (np.log2(1 + snr_main) - 
                 np.log2(1 + snr_best_eve))
            
            if cs < target_rate:
                outage_count += 1
        
        return outage_count / n_trials
    
    def visualize_network(self, area_radius=100):
        """Visualize a random network deployment"""
        network = self.deploy_network(area_radius)
        
        plt.figure(figsize=(8, 8))
        
        # Plot eavesdroppers
        if network['n_eves'] > 0:
            plt.scatter(
                network['eve_positions'][:, 0],
                network['eve_positions'][:, 1],
                c='red', marker='x', s=100,
                label=f'Eavesdroppers (n={network["n_eves"]})',
                zorder=3)
        
        # Plot transmitter
        plt.scatter([0], [0], c='blue', marker='^', 
                   s=200, label='Transmitter', zorder=4)
        
        # Plot legitimate receiver
        plt.scatter([10], [0], c='green', marker='o',
                   s=200, label='Legitimate Receiver', 
                   zorder=4)
        
        circle = plt.Circle((0, 0), area_radius,
                            fill=False, color='gray',
                            linestyle='--', alpha=0.5)
        plt.gca().add_patch(circle)
        
        plt.xlim(-area_radius*1.1, area_radius*1.1)
        plt.ylim(-area_radius*1.1, area_radius*1.1)
        plt.legend()
        plt.title(f'Random Network — λ={self.lambda_e}/m²')
        plt.grid(True, alpha=0.3)
        plt.axis('equal')
        plt.show()


# Analysis
model = StochasticGeometryPLS(
    lambda_eve=0.001,
    tx_power=1.0,
    noise_power=0.01,
    path_loss_exp=3.5
)

# Visualize network deployment
model.visualize_network()

# SOP vs eavesdropper density
densities = np.logspace(-4, -2, 10)
sop_values = []

for density in densities:
    m = StochasticGeometryPLS(density, 1.0, 0.01, 3.5)
    sop = m.monte_carlo_sop(
        target_rate=1.0, n_trials=500)
    sop_values.append(sop)
    print(f"λ={density:.4f}: SOP={sop:.3f}")

plt.figure(figsize=(10, 5))
plt.semilogx(densities, sop_values, 
             'o-', linewidth=2)
plt.xlabel('Eavesdropper Density (per m²)')
plt.ylabel('Secrecy Outage Probability')
plt.title('SOP vs Eavesdropper Density — PPP Model')
plt.grid(True, alpha=0.3)
plt.show()

What the code does: This class models a realistic wireless network where eavesdroppers are scattered randomly according to a Poisson Point Process. The deploy_network() method generates a random realization — sometimes there are 3 eavesdroppers nearby, sometimes 15, sometimes none. The monte_carlo_sop() method runs hundreds of these random deployments and computes how often secrecy fails across all of them. The visualization shows you what one random network snapshot looks like — the transmitter in blue, the legitimate receiver in green, and randomly scattered eavesdroppers in red.

What the math means: In real wireless networks, we do not know where eavesdroppers are. They could be anywhere. Stochastic geometry models this uncertainty by treating eavesdropper locations as a random spatial process. The Poisson Point Process is the natural choice — it models completely random spatial distributions with no clustering or repulsion. The density parameter controls how many eavesdroppers exist on average per unit area. As density increases, the probability that at least one eavesdropper has a strong channel to the transmitter also increases, driving up the secrecy outage probability.

Connection to Secrecy Analysis:

For a PPP of eavesdroppers with density $\lambda_e$λe​, the secrecy outage probability can be expressed as:$$P_{out} = 1 – \mathcal{L}_{\Phi_e}\left(\frac{\phi}{\bar{\gamma}_m}\right)$$Where $\mathcal{L}_{\Phi_e}(\cdot)$ is the Laplace functional of the eavesdropper point process and $\phi = 2^{R_s} – 1$ is the secrecy rate threshold.

---

Probabilistic Models + Kalman Filter

python

import numpy as np
import matplotlib.pyplot as plt

class KalmanChannelEstimator:
    """
    Kalman filter for secure wireless channel estimation
    Tracks channel state for Physical Layer Security
    """
    
    def __init__(self, process_noise=0.01, 
                 measurement_noise=0.1,
                 initial_state=1.0,
                 initial_uncertainty=1.0):
        """
        Args:
            process_noise: Channel variation variance (Q)
            measurement_noise: Observation noise variance (R)
            initial_state: Initial channel estimate
            initial_uncertainty: Initial estimation uncertainty
        """
        self.Q = process_noise
        self.R = measurement_noise
        self.x = initial_state
        self.P = initial_uncertainty
        
        # History storage
        self.estimates = [initial_state]
        self.uncertainties = [initial_uncertainty]
        self.gains = []
    
    def predict(self):
        """
        Prediction step — propagate state forward
        Channel is assumed to evolve slowly (AR model)
        """
        # State prediction (channel stays similar)
        x_pred = self.x
        
        # Covariance prediction (uncertainty grows)
        P_pred = self.P + self.Q
        
        return x_pred, P_pred
    
    def update(self, measurement):
        """
        Update step — incorporate new measurement
        
        Args:
            measurement: Noisy channel observation
        
        Returns:
            Updated state estimate and uncertainty
        """
        # Predict
        x_pred, P_pred = self.predict()
        
        # Kalman gain
        K = P_pred / (P_pred + self.R)
        
        # State update
        self.x = x_pred + K * (measurement - x_pred)
        
        # Covariance update
        self.P = (1 - K) * P_pred
        
        # Store history
        self.estimates.append(self.x)
        self.uncertainties.append(self.P)
        self.gains.append(K)
        
        return self.x, self.P, K
    
    def estimate_secrecy_capacity(self, 
                                   eve_estimate,
                                   tx_snr=10.0):
        """
        Estimate secrecy capacity from channel estimates
        
        Args:
            eve_estimate: Estimated eavesdropper channel
            tx_snr: Transmit SNR
        
        Returns:
            Estimated secrecy capacity
        """
        snr_main = tx_snr * self.x**2
        snr_eve = tx_snr * eve_estimate**2
        return max(0, np.log2(1 + snr_main) - 
                      np.log2(1 + snr_eve))


# Simulate channel tracking
np.random.seed(42)
n_steps = 200

# True channel evolution (slow fading)
true_channel = np.zeros(n_steps)
true_channel[0] = 1.0
for t in range(1, n_steps):
    true_channel[t] = (0.99 * true_channel[t-1] + 
                       np.random.normal(0, 0.05))

# Noisy measurements
measurements = true_channel + np.random.normal(
    0, 0.3, n_steps)

# Kalman filter estimation
kf_main = KalmanChannelEstimator(
    process_noise=0.01,
    measurement_noise=0.09,
    initial_state=measurements[0]
)

# Eavesdropper channel estimation
true_eve_channel = np.zeros(n_steps)
true_eve_channel[0] = 0.7
for t in range(1, n_steps):
    true_eve_channel[t] = (0.99 * true_eve_channel[t-1] + 
                           np.random.normal(0, 0.04))

eve_measurements = true_eve_channel + np.random.normal(
    0, 0.3, n_steps)

kf_eve = KalmanChannelEstimator(
    process_noise=0.01,
    measurement_noise=0.09,
    initial_state=eve_measurements[0]
)

# Run filters
secrecy_estimates = []
for t in range(1, n_steps):
    main_est, main_unc, main_K = kf_main.update(
        measurements[t])
    eve_est, eve_unc, eve_K = kf_eve.update(
        eve_measurements[t])
    
    cs_est = kf_main.estimate_secrecy_capacity(
        eve_est, tx_snr=10.0)
    secrecy_estimates.append(cs_est)

# True secrecy capacity
true_cs = np.maximum(0, 
    np.log2(1 + 10*true_channel**2) - 
    np.log2(1 + 10*true_eve_channel**2))

# Plotting
fig, axes = plt.subplots(3, 1, figsize=(12, 12))

# Channel estimates
axes[0].plot(true_channel, label='True Main Channel',
             linewidth=1.5, alpha=0.8)
axes[0].plot(kf_main.estimates[:-1], '--',
             label='Kalman Estimate', linewidth=1.5)
axes[0].plot(measurements, '.', alpha=0.3,
             label='Noisy Measurements', markersize=2)
axes[0].set_ylabel('Channel Magnitude')
axes[0].set_title('Kalman Filter Channel Tracking')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# Kalman gain evolution
axes[1].plot(kf_main.gains, linewidth=1.5,
             color='orange')
axes[1].set_ylabel('Kalman Gain')
axes[1].set_title('Kalman Gain Evolution')
axes[1].grid(True, alpha=0.3)

# Secrecy capacity estimation
axes[2].plot(true_cs[1:], label='True Secrecy Capacity',
             linewidth=1.5, alpha=0.8)
axes[2].plot(secrecy_estimates, '--',
             label='Estimated Secrecy Capacity',
             linewidth=1.5)
axes[2].set_ylabel('Secrecy Capacity (bits/use)')
axes[2].set_xlabel('Time Step')
axes[2].set_title('Secrecy Capacity Estimation via Kalman Filter')
axes[2].legend()
axes[2].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

print(f"Channel estimation MSE: "
      f"{np.mean((true_channel[1:] - kf_main.estimates[:-2])**2):.6f}")
print(f"Secrecy capacity MSE: "
      f"{np.mean((true_cs[1:] - secrecy_estimates)**2):.6f}")

What the code does: This implementation tracks two wireless channels simultaneously using separate Kalman filters — one for the legitimate channel and one for the estimated eavesdropper channel. At each time step, the filter takes a noisy measurement of the current channel state, runs the predict-update cycle, and produces an improved estimate. The secrecy capacity is then computed from these estimates rather than from the raw noisy measurements. The three plots show the channel tracking quality, how the Kalman gain stabilizes over time, and how accurately the secrecy capacity can be estimated.

What the math means: In practice, the transmitter never knows the exact channel state. It can only measure a noisy version of it. Without filtering, using raw measurements to make security decisions leads to poor performance — sometimes thinking the channel is secure when it is not, sometimes thinking it is insecure when it actually is. The Kalman filter solves this by maintaining a running estimate of the channel state that is statistically optimal — no algorithm can do better given the noise statistics. The gain K automatically balances trust in the prediction versus trust in the new measurement, just as we described in the Kalman filter article.

State Space Model for Channel Tracking:$$h[t+1] = \alpha h[t] + w[t], \quad w[t] \sim \mathcal{N}(0, Q)$$$$y[t] = h[t] + v[t], \quad v[t] \sim \mathcal{N}(0, R)$$Where $\alpha$ is the channel correlation coefficient capturing how slowly the channel changes, $Q$ is the process noise variance controlling channel variability, and $R$ is the measurement noise variance from pilot-based estimation.

---

Machine Learning + Probabilistic PLS Models

python

import numpy as np
import matplotlib.pyplot as plt
from sklearn.neural_network import MLPRegressor
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import train_test_split

class DeepLearningPLS:
    """
    Neural network for PLS secrecy rate prediction
    and optimization
    """
    
    def __init__(self, hidden_layers=(128, 64, 32)):
        """
        Args:
            hidden_layers: Tuple defining NN architecture
        """
        self.model = MLPRegressor(
            hidden_layer_sizes=hidden_layers,
            activation='relu',
            max_iter=500,
            random_state=42
        )
        self.scaler_X = StandardScaler()
        self.scaler_y = StandardScaler()
        self.is_trained = False
    
    def generate_training_data(self, n_samples=10000):
        """
        Generate synthetic PLS training data
        
        Features: channel statistics, system parameters
        Target: optimal secrecy rate
        """
        # Random system parameters
        mean_snr_main = np.random.uniform(1, 20, n_samples)
        mean_snr_eve = np.random.uniform(0.5, 10, n_samples)
        path_loss_exp = np.random.uniform(2, 4, n_samples)
        n_antennas = np.random.randint(1, 8, n_samples)
        lambda_eve = np.random.uniform(0.0001, 0.01, n_samples)
        
        features = np.column_stack([
            mean_snr_main,
            mean_snr_eve,
            path_loss_exp,
            n_antennas,
            lambda_eve,
            mean_snr_main / mean_snr_eve,  # SNR ratio
            np.log(mean_snr_main),
            np.log(mean_snr_eve)
        ])
        
        # Compute optimal secrecy rates
        optimal_rates = []
        for i in range(n_samples):
            # Monte Carlo for each configuration
            h_m = np.random.rayleigh(
                np.sqrt(mean_snr_main[i]/2), 1000)
            h_e = np.random.rayleigh(
                np.sqrt(mean_snr_eve[i]/2), 1000)
            
            snr_m = h_m**2
            snr_e = h_e**2
            
            # Find optimal rate
            best_throughput = 0
            best_rate = 0
            for r in np.linspace(0.1, 5.0, 50):
                cs = np.maximum(0, 
                    np.log2(1 + snr_m) - 
                    np.log2(1 + snr_e))
                throughput = r * np.mean(cs >= r)
                if throughput > best_throughput:
                    best_throughput = throughput
                    best_rate = r
            
            optimal_rates.append(best_rate)
        
        return features, np.array(optimal_rates)
    
    def train(self, n_samples=5000):
        """Train the neural network"""
        print("Generating training data...")
        X, y = self.generate_training_data(n_samples)
        
        X_train, X_test, y_train, y_test = \
            train_test_split(X, y, test_size=0.2)
        
        X_train_scaled = self.scaler_X.fit_transform(X_train)
        X_test_scaled = self.scaler_X.transform(X_test)
        
        print("Training neural network...")
        self.model.fit(X_train_scaled, y_train)
        self.is_trained = True
        
        train_score = self.model.score(
            X_train_scaled, y_train)
        test_score = self.model.score(
            X_test_scaled, y_test)
        
        print(f"Train R² score: {train_score:.4f}")
        print(f"Test R²  score: {test_score:.4f}")
        
        return test_score
    
    def predict_optimal_rate(self, mean_snr_main, 
                              mean_snr_eve,
                              path_loss_exp=3.5,
                              n_antennas=1,
                              lambda_eve=0.001):
        """
        Predict optimal secrecy rate for given parameters
        """
        if not self.is_trained:
            raise ValueError("Model not trained yet")
        
        features = np.array([[
            mean_snr_main,
            mean_snr_eve,
            path_loss_exp,
            n_antennas,
            lambda_eve,
            mean_snr_main / mean_snr_eve,
            np.log(mean_snr_main),
            np.log(mean_snr_eve)
        ]])
        
        features_scaled = self.scaler_X.transform(features)
        return max(0, self.model.predict(features_scaled)[0])


# Train and test
dl_pls = DeepLearningPLS(hidden_layers=(64, 32, 16))
dl_pls.train(n_samples=3000)

# Test predictions
test_cases = [
    (10, 3, "Strong main, weak eve"),
    (5,  5, "Equal channels"),
    (3, 10, "Weak main, strong eve"),
    (15, 2, "Very strong advantage"),
]

print("\nOptimal Secrecy Rate Predictions:")
print("-" * 50)
for snr_m, snr_e, desc in test_cases:
    rate = dl_pls.predict_optimal_rate(snr_m, snr_e)
    print(f"{desc}: {rate:.3f} bits/channel use")

What the code does: This neural network learns to predict the optimal secrecy transmission rate given system parameters, without needing to run expensive Monte Carlo simulations at deployment time. The training phase generates thousands of random system configurations, computes the optimal rate for each using Monte Carlo simulation, and trains the network to map from system parameters to optimal rate. Once trained, predictions take microseconds instead of seconds. This is practically valuable in dynamic wireless systems where channel conditions change frequently and real-time rate adaptation is needed.

What the math means: Traditional PLS optimization requires solving integrals or running simulations for every new channel configuration. In a dynamic network, this is too slow. Deep learning amortizes this computation — the expensive optimization happens once during training, and the neural network serves as a fast approximator. This is the machine learning approach to probabilistic optimization: replace an expensive computation with a learned function that approximates it. The key insight is that the mapping from channel statistics to optimal security parameters, while complex, has structure that a neural network can learn and generalize.$$\hat{R}_s^* = f_\theta(\bar{\gamma}_m, \bar{\gamma}_e, \alpha, N_t, \lambda_e)$$Where $f_\theta$​ is the neural network with parameters $\theta$θ learned to approximate the true optimal rate $R_s^* = \arg\max_{R_s} \tau_s(R_s)$.

---

Probabilistic Control Engineering for PLS

python

import numpy as np
import matplotlib.pyplot as plt

class PCESecurityController:
    """
    Probabilistic Control Engineering framework
    applied to Physical Layer Security
    
    Treats secure transmission as a control problem:
    - State: current channel quality difference
    - Control: transmission rate and power allocation
    - Objective: maintain secrecy while maximizing throughput
    """
    
    def __init__(self, target_sop=0.05, 
                 learning_rate=0.01):
        """
        Args:
            target_sop: Target secrecy outage probability
            learning_rate: Rate adaptation step size
        """
        self.target_sop = target_sop
        self.lr = learning_rate
        
        # Controller state
        self.current_rate = 1.0
        self.estimated_sop = 0.5
        
        # History
        self.rate_history = [self.current_rate]
        self.sop_history = [self.estimated_sop]
        self.cs_history = []
        
        # Kalman filter for channel tracking
        self.channel_estimate = 1.0
        self.channel_uncertainty = 1.0
        self.Q = 0.01  # Process noise
        self.R = 0.09  # Measurement noise
    
    def update_channel_estimate(self, measurement):
        """Kalman filter update for channel state"""
        # Predict
        P_pred = self.channel_uncertainty + self.Q
        
        # Update
        K = P_pred / (P_pred + self.R)
        self.channel_estimate = (self.channel_estimate + 
            K * (measurement - self.channel_estimate))
        self.channel_uncertainty = (1 - K) * P_pred
        
        return self.channel_estimate
    
    def estimate_sop(self, snr_main, snr_eve, 
                      window=50):
        """
        Online SOP estimation from recent observations
        """
        cs = max(0, np.log2(1 + snr_main) - 
                    np.log2(1 + snr_eve))
        self.cs_history.append(cs)
        
        if len(self.cs_history) > window:
            recent_cs = self.cs_history[-window:]
            return np.mean(
                np.array(recent_cs) < self.current_rate)
        return self.estimated_sop
    
    def adapt_rate(self):
        """
        Control law: adapt rate based on SOP feedback
        
        If SOP > target: decrease rate (safer)
        If SOP < target: increase rate (more efficient)
        """
        error = self.estimated_sop - self.target_sop
        
        # Proportional control
        rate_adjustment = -self.lr * error
        
        # Update rate with constraints
        self.current_rate = np.clip(
            self.current_rate + rate_adjustment,
            0.1, 5.0)
        
        return self.current_rate
    
    def run_simulation(self, n_steps=500, 
                        mean_snr_main=10,
                        mean_snr_eve=3):
        """
        Run closed-loop PLS control simulation
        """
        for t in range(n_steps):
            # Channel realizations
            h_main = np.random.rayleigh(
                np.sqrt(mean_snr_main/2))
            h_eve = np.random.rayleigh(
                np.sqrt(mean_snr_eve/2))
            
            snr_main = h_main**2
            snr_eve = h_eve**2
            
            # Update channel estimate
            self.update_channel_estimate(snr_main)
            
            # Estimate current SOP
            self.estimated_sop = self.estimate_sop(
                snr_main, snr_eve)
            
            # Adapt transmission rate
            self.current_rate = self.adapt_rate()
            
            # Store history
            self.rate_history.append(self.current_rate)
            self.sop_history.append(self.estimated_sop)
        
        return self.rate_history, self.sop_history


# Run PCE controller
controller = PCESecurityController(
    target_sop=0.05,
    learning_rate=0.05
)

rates, sops = controller.run_simulation(
    n_steps=500, 
    mean_snr_main=10,
    mean_snr_eve=3
)

fig, axes = plt.subplots(2, 1, figsize=(12, 8))

axes[0].plot(rates, linewidth=1.5, color='steelblue',
             label='Adapted Rate')
axes[0].axhline(y=1.0, color='gray', linestyle='--',
                label='Initial Rate', alpha=0.5)
axes[0].set_ylabel('Secrecy Rate (bits/use)')
axes[0].set_title('PCE Rate Adaptation for PLS')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

axes[1].plot(sops, linewidth=1.5, color='orange',
             label='Estimated SOP')
axes[1].axhline(y=0.05, color='red', linestyle='--',
                label='Target SOP = 0.05', linewidth=2)
axes[1].set_ylabel('Secrecy Outage Probability')
axes[1].set_xlabel('Time Step')
axes[1].set_title('SOP Convergence to Target')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

final_sop = np.mean(sops[-100:])
final_rate = np.mean(rates[-100:])
print(f"Final average SOP: {final_sop:.4f} "
      f"(target: 0.05)")
print(f"Final average rate: {final_rate:.4f} bits/use")

What the code does: This controller treats Physical Layer Security as a feedback control problem — exactly the Probabilistic Control Engineering framework from our earlier articles. The system continuously measures how often secrecy outages occur, compares that to the target outage probability, computes the error, and adjusts the transmission rate accordingly. If too many outages occur, the rate drops to be more conservative. If outages are rare, the rate increases to capture more throughput. The Kalman filter runs inside the controller to maintain a clean channel estimate from noisy measurements. Watch the plots — the SOP starts far from the target and gradually converges as the controller learns the right rate.

What the math means: This is the PCE framework applied to wireless security. The channel state is random and uncertain — handled by the Kalman filter. The security objective is probabilistic — we target a specific SOP rather than a specific channel condition. The control action is the transmission rate — adjusted continuously based on feedback from observed security performance. This closed-loop approach automatically adapts to changing channel conditions, eavesdropper behavior, and network dynamics without requiring prior knowledge of the exact statistical distributions. It is robust in the same way that feedback control is always more robust than open-loop control.

PCE Control Law for PLS:$$R_s[t+1] = R_s[t] – \mu \cdot (P_{out}[t] – P_{out}^{target})$$Where $\mu$ is the adaptation step size, $P_{out}[t]$ is the estimated secrecy outage probability at time $t$, and $P_{out}^{target}$ is the desired security level. This proportional controller drives the system toward the target SOP in expectation.

---

RIS — Reconfigurable Intelligent Surfaces

python

import numpy as np
import matplotlib.pyplot as plt
from itertools import product

class RISAssistedPLS:
    """
    RIS-Assisted Physical Layer Security
    Optimizes RIS phase shifts for secrecy maximization
    """
    
    def __init__(self, n_elements=16, 
                 tx_power=1.0,
                 noise_power=0.01):
        """
        Args:
            n_elements: Number of RIS reflecting elements
            tx_power: Transmitter power
            noise_power: Noise power at receivers
        """
        self.N = n_elements
        self.P = tx_power
        self.N0 = noise_power
        
        # Random phase shifts initially
        self.phase_shifts = np.random.uniform(
            0, 2*np.pi, n_elements)
    
    def generate_channels(self, 
                          d_tx_ris=20.0,
                          d_ris_rx=15.0,
                          d_ris_eve=25.0,
                          d_direct_rx=50.0,
                          d_direct_eve=45.0):
        """
        Generate channel coefficients for all links
        
        Returns:
            Dictionary of complex channel vectors
        """
        path_loss = lambda d: d**(-3.5)
        
        # TX-RIS channel (N x 1 vector)
        h_tx_ris = (np.sqrt(path_loss(d_tx_ris)/2) * 
                   (np.random.randn(self.N) + 
                    1j*np.random.randn(self.N)))
        
        # RIS-RX channel (1 x N vector)
        h_ris_rx = (np.sqrt(path_loss(d_ris_rx)/2) * 
                   (np.random.randn(self.N) + 
                    1j*np.random.randn(self.N)))
        
        # RIS-Eve channel (1 x N vector)
        h_ris_eve = (np.sqrt(path_loss(d_ris_eve)/2) * 
                    (np.random.randn(self.N) + 
                     1j*np.random.randn(self.N)))
        
        # Direct TX-RX channel
        h_direct_rx = (np.sqrt(path_loss(d_direct_rx)/2) * 
                      (np.random.randn() + 
                       1j*np.random.randn()))
        
        # Direct TX-Eve channel
        h_direct_eve = (np.sqrt(path_loss(d_direct_eve)/2) * 
                       (np.random.randn() + 
                        1j*np.random.randn()))
        
        return {
            'h_tx_ris': h_tx_ris,
            'h_ris_rx': h_ris_rx,
            'h_ris_eve': h_ris_eve,
            'h_direct_rx': h_direct_rx,
            'h_direct_eve': h_direct_eve
        }
    
    def compute_effective_snr(self, channels):
        """
        Compute SNR at legitimate receiver and eavesdropper
        considering both direct and RIS-reflected paths
        """
        # RIS phase shift matrix
        Phi = np.diag(np.exp(1j * self.phase_shifts))
        
        # Effective channel gains
        # RIS path: h_ris_rx * Phi * h_tx_ris
        ris_contrib_rx = (channels['h_ris_rx'] @ 
                         Phi @ channels['h_tx_ris'])
        
        ris_contrib_eve = (channels['h_ris_eve'] @ 
                          Phi @ channels['h_tx_ris'])
        
        # Total effective channels
        h_eff_rx = channels['h_direct_rx'] + ris_contrib_rx
        h_eff_eve = channels['h_direct_eve'] + ris_contrib_eve
        
        # SNRs
        snr_rx = (self.P * np.abs(h_eff_rx)**2 / self.N0)
        snr_eve = (self.P * np.abs(h_eff_eve)**2 / self.N0)
        
        return snr_rx, snr_eve
    
    def optimize_phases_random(self, 
                                n_trials=1000,
                                n_channel_samples=100):
        """
        Random search phase optimization for secrecy
        """
        best_cs = -np.inf
        best_phases = self.phase_shifts.copy()
        
        for trial in range(n_trials):
            # Random phase configuration
            test_phases = np.random.uniform(
                0, 2*np.pi, self.N)
            self.phase_shifts = test_phases
            
            # Average secrecy over channel realizations
            cs_samples = []
            for _ in range(n_channel_samples):
                channels = self.generate_channels()
                snr_rx, snr_eve = self.compute_effective_snr(
                    channels)
                cs = max(0, np.log2(1 + snr_rx) - 
                           np.log2(1 + snr_eve))
                cs_samples.append(cs)
            
            avg_cs = np.mean(cs_samples)
            
            if avg_cs > best_cs:
                best_cs = avg_cs
                best_phases = test_phases.copy()
        
        self.phase_shifts = best_phases
        return best_phases, best_cs
    
    def compare_with_without_ris(self, n_samples=1000):
        """
        Compare secrecy performance with and without RIS
        """
        # Without RIS
        cs_no_ris = []
        cs_with_ris = []
        
        for _ in range(n_samples):
            channels = self.generate_channels()
            
            # No RIS — direct path only
            snr_rx_direct = (self.P * 
                np.abs(channels['h_direct_rx'])**2 / 
                self.N0)
            snr_eve_direct = (self.P * 
                np.abs(channels['h_direct_eve'])**2 / 
                self.N0)
            cs_no_ris.append(max(0, 
                np.log2(1 + snr_rx_direct) - 
                np.log2(1 + snr_eve_direct)))
            
            # With RIS
            snr_rx, snr_eve = self.compute_effective_snr(
                channels)
            cs_with_ris.append(max(0, 
                np.log2(1 + snr_rx) - 
                np.log2(1 + snr_eve)))
        
        return np.array(cs_no_ris), np.array(cs_with_ris)


# RIS Analysis
ris_system = RISAssistedPLS(n_elements=16)

print("Optimizing RIS phase shifts...")
best_phases, best_cs = ris_system.optimize_phases_random(
    n_trials=200, n_channel_samples=50)
print(f"Optimized Average Secrecy Rate: {best_cs:.4f} bits/use")

print("\nComparing with and without RIS...")
cs_no_ris, cs_with_ris = ris_system.compare_with_without_ris(
    n_samples=500)

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# CDF comparison
for cs_arr, label, color in [
    (cs_no_ris, 'Without RIS', 'red'),
    (cs_with_ris, 'With Optimized RIS', 'steelblue')]:
    cs_sorted = np.sort(cs_arr)
    cdf = np.arange(1, len(cs_sorted)+1) / len(cs_sorted)
    axes[0].plot(cs_sorted, cdf, linewidth=2,
                label=label, color=color)

axes[0].set_xlabel('Secrecy Capacity (bits/use)')
axes[0].set_ylabel('CDF')
axes[0].set_title('RIS Impact on Secrecy Capacity CDF')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# Bar comparison
metrics = ['Mean Cs', 'P(Cs>0)', 'P(Cs>1)']
no_ris_vals = [
    np.mean(cs_no_ris),
    np.mean(cs_no_ris > 0),
    np.mean(cs_no_ris > 1)
]
ris_vals = [
    np.mean(cs_with_ris),
    np.mean(cs_with_ris > 0),
    np.mean(cs_with_ris > 1)
]

x = np.arange(len(metrics))
width = 0.35
axes[1].bar(x - width/2, no_ris_vals, width,
           label='Without RIS', color='red', alpha=0.7)
axes[1].bar(x + width/2, ris_vals, width,
           label='With RIS', color='steelblue', alpha=0.7)
axes[1].set_xticks(x)
axes[1].set_xticklabels(metrics)
axes[1].set_ylabel('Value')
axes[1].set_title('PLS Metrics: With vs Without RIS')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

print(f"\nWithout RIS — Mean Cs: {np.mean(cs_no_ris):.4f}")
print(f"With RIS    — Mean Cs: {np.mean(cs_with_ris):.4f}")
print(f"Improvement: {(np.mean(cs_with_ris)/max(np.mean(cs_no_ris),0.001)-1)*100:.1f}%")

What the code does: This implementation models a Reconfigurable Intelligent Surface — a programmable array of reflecting elements that can be tuned to redirect wireless signals. The optimize_phases_random() method tries hundreds of random phase configurations and keeps the one that produces the highest average secrecy capacity. The compare_with_without_ris() method then shows clearly how much the optimized RIS improves security. The CDF plot shows the entire distribution shift — with a good RIS configuration, not just the average but the entire secrecy capacity distribution moves upward.

What the math means: A Reconfigurable Intelligent Surface works by reflecting incoming wireless signals with controlled phase shifts. By carefully tuning these phase shifts, the reflected signal from the RIS can be made to add constructively at the legitimate receiver and destructively at the eavesdropper. This is beamforming extended to an intelligent reflecting surface. The key probabilistic insight is that RIS optimization changes the underlying distribution of the effective channel — not just a single channel realization. A well-optimized RIS systematically shifts the entire distribution of the legitimate channel upward while keeping the eavesdropper’s channel distribution unchanged.

RIS-Assisted Effective Channel:$$h_{eff} = h_{direct} + \mathbf{h}_{RIS-RX}^T \boldsymbol{\Phi} \mathbf{h}_{TX-RIS}$$Where $\boldsymbol{\Phi} = \text{diag}(e^{j\theta_1}, \ldots, e^{j\theta_N})$ is the RIS phase shift matrix with $\theta_n \in [0, 2\pi)$ for each of the $N$ reflecting elements.

---

Future of PLS — 6G and Beyond

The trajectory of Physical Layer Security research points clearly toward three converging developments in the coming years.

AI-native PLS systems will become standard. Rather than designing fixed probabilistic models and then adding machine learning on top, future systems will be designed from the ground up with learned probabilistic representations. The channel model, the security metric, and the optimization algorithm will all be jointly learned from data rather than derived analytically.

Terahertz band communications, central to 6G systems, will create new PLS opportunities and challenges. The extremely high frequencies and short wavelengths of THz bands produce fading statistics that differ significantly from sub-6 GHz models. Channel sparsity in THz bands — where only a few dominant paths exist — creates strong opportunities for PLS because the legitimate and eavesdropper channels are even more physically distinct. But it also creates new modeling challenges that the Rayleigh and Nakagami frameworks do not fully capture.

Quantum Physical Layer Security represents the ultimate convergence — using quantum mechanical properties of electromagnetic fields to achieve security guarantees that go beyond classical information theory. Quantum key distribution at the physical layer is no longer purely theoretical, and within a decade it may become practical for short-range wireless links.

---

Conclusion

Physical Layer Security is one of those fields where classical engineering theory, modern probabilistic modeling, and cutting-edge machine learning all converge in a genuinely meaningful way.

The secrecy capacity framework from information theory tells us what is theoretically possible. Probabilistic fading models — Rayleigh, Nakagami, Rician — tell us how channel randomness affects achievable security. Stochastic geometry tells us how the spatial distribution of eavesdroppers shapes network-level security. Kalman filtering bridges the classical control theory tradition to real-time channel tracking for security decisions. And deep learning is beginning to automate the optimization problems that were previously too complex for analytical solutions.

The PCE framework ties all of this together. Treating Physical Layer Security as a probabilistic control problem — with uncertain channel states, stochastic objectives, and adaptive feedback control — gives us both the theoretical foundation and the practical tools to build wireless systems that are secure by physics, not just by computational hardness.

For AI cybersecurity engineers, PLS represents a frontier where your skills in probabilistic modeling, machine learning, and systems thinking are genuinely needed. The mathematics is rich, the applications are real, and the problems are open.