如何检测NumPy数组中SIS模型模拟的未知双值稳态波动？

Great question! Detecting that persistent, bounded fluctuation in your SIS model time series is crucial for identifying when the infection has reached an endemic steady state (instead of dying out or still growing). Here’s a practical, code-driven approach tailored to your scenario:

First, Clarify Key Thresholds #

Let’s start by formalizing the thresholds you mentioned:

• Extinction Threshold: Infection count drops below this value → dies out.
• Trigger Threshold: 40% coverage (infected nodes / total nodes) → start checking for steady-state fluctuations.
• Steady State: Infection bounces between two unknown values without trending up/down.

Step-by-Step Detection Implementation #

We’ll build on your sample simulation code and add detection logic using rolling statistics, stationarity tests, and even clustering for precise bounds.

1. Full Code with Detection #

import random
import numpy as np
import matplotlib.pyplot as plt
from statsmodels.tsa.stattools import adfuller
from sklearn.cluster import KMeans

# Simulation parameters
total_nodes = 50  # Total nodes in your network
initial_steps = 10
fluct_min = 11
fluct_max = 17
extinction_threshold = 5  # Infection count below this = die out
coverage_trigger = 0.4  # 40% coverage to start checking

# Generate simulated SIS-like time series (matches your sample)
infection_counts = np.arange(0, initial_steps, 1)
for _ in range(100):
    infection_counts = np.append(infection_counts, random.randint(fluct_min, fluct_max))

# Calculate coverage over time
coverage = infection_counts / total_nodes

2. Isolate the Relevant Time Window #

First, filter out transient phases and focus on data where coverage crosses your 40% trigger:

# Find the first time coverage hits 40% and stays above extinction threshold
trigger_indices = np.where((coverage >= coverage_trigger) & (infection_counts >= extinction_threshold))[0]

if len(trigger_indices) == 0:
    print("Infection never reached the 40% coverage threshold, or died out early.")
else:
    # Focus on the steady-state candidate segment
    steady_candidate = infection_counts[trigger_indices[0]:]
    steady_coverage = coverage[trigger_indices[0]:]

3. Detect Steady-State Fluctuations #

We’ll use three complementary checks to confirm the fluctuation phase:

A. Rolling Statistics Check

This verifies that the infection count isn’t trending—its min/max values stabilize over time:

window_size = 10  # Adjust based on your time step frequency

# Calculate rolling min/max to spot stable bounds
rolling_min = np.convolve(steady_candidate, np.ones(window_size)/window_size, mode='valid')
rolling_max = np.convolve(steady_candidate, np.ones(window_size)/window_size, mode='valid')

# Check if rolling stats have low variance (no upward/downward trend)
min_variance = np.var(rolling_min)
max_variance = np.var(rolling_max)
is_stable_bounds = (min_variance < 1) and (max_variance < 1)  # Tune variance threshold to your data

B. Stationarity Test (ADF)

The Augmented Dickey-Fuller test checks if the time series is stationary (mean/variance don’t change over time)—a hallmark of steady-state fluctuation:

adf_result = adfuller(steady_candidate)
p_value = adf_result[1]
is_stationary = p_value < 0.05  # Low p-value = stationary

C. Cluster Detection for Fluctuation Bounds

If your steady state strictly bounces between two values, K-means clustering will pinpoint those exact bounds:

# Reshape data for clustering
cluster_data = steady_candidate.reshape(-1, 1)
kmeans = KMeans(n_clusters=2, random_state=0).fit(cluster_data)
fluct_bounds = sorted(kmeans.cluster_centers_.flatten())

4. Combine Checks and Visualize #

# Final decision logic
if is_stable_bounds and is_stationary:
    print("✅ Infection has entered a steady-state fluctuation phase!")
    print(f"Fluctuating between ~{fluct_bounds[0]:.2f} and {fluct_bounds[1]:.2f} infected nodes.")
else:
    print("❌ Infection is not in a steady-state fluctuation (still growing/decaying).")

# Plot results for verification
plt.figure(figsize=(10, 6))
plt.plot(infection_counts, label='Infected Nodes')
plt.axhline(y=extinction_threshold, color='r', linestyle='--', label='Extinction Threshold')
plt.axhline(y=coverage_trigger * total_nodes, color='g', linestyle='--', label='40% Coverage Trigger')
plt.axvline(x=trigger_indices[0], color='purple', linestyle=':', label='Steady-State Start')
plt.xlabel('Time Step')
plt.ylabel('Number of Infected Nodes')
plt.legend()
plt.title('SIS Model Infection Count with Steady-State Detection')
plt.show()

Key Tips for Tuning #

• Transient Phase: Add a check to ensure the series stays above the trigger threshold for multiple consecutive steps (e.g., 5 time steps) to avoid false positives from temporary spikes.
• Window Size: Adjust the rolling window size based on your simulation’s time step frequency—larger windows smooth noise but might miss rapid fluctuations.
• Variance Threshold: Tweak the variance cutoff in the rolling stats check to match your model’s noise level (higher infection/recovery noise may require a higher threshold).

内容的提问来源于stack exchange，提问作者Terry BRETT