import numpy as np
 import scipy.stats as stats
 import matplotlib.pyplot as plt

# Simulating normal system behavior (3 correlated sensors)
np.random.seed(42)
mean = [50, 75, 100]  # Mean values for three sensors
cov = [[10, 5, 2], [5, 15, 3], [2, 3, 20]]  # Covariance matrix

# Generate 500 normal operation samples
 normal_data = np.random.multivariate_normal(mean, cov, size=500)

class MSET:
    def __init__(self, memory_matrix):
        self.memory_matrix = memory_matrix  # Store normal system states

    def estimate(self, input_vector):
        """
        Estimates the expected state based on historical data.
        Uses nearest neighbors to compute weighted estimation.
        """
        weights = np.exp(-np.linalg.norm(self.memory_matrix - input_vector, axis=1))
        weights /= np.sum(weights)
         return np.dot(weights, self.memory_matrix)

# Initialize MSET with normal data as memory
 mset_model = MSET(memory_matrix=normal_data)

# filepath: deephub\5\20250327\article.md
# Simulated test data (normal + some anomalies)
test_data = np.vstack([
    np.random.multivariate_normal(mean, cov, size=450),  # Normal
    np.random.multivariate_normal([70, 50, 130], cov, size=50)  # Anomalies
])

# Compute estimated values
estimated_data = np.array([mset_model.estimate(x) for x in test_data])
# Compute residuals
 residuals = np.linalg.norm(test_data - estimated_data, axis=1)

# Define thresholds for SPRT
alpha = 0.05  # False positive rate
beta = 0.05   # False negative rate
mu_0, sigma_0 = np.mean(residuals[:450]), np.std(residuals[:450])  # Normal behavior
mu_1 = mu_0 + 3 * sigma_0  # Anomalous mean shift

# SPRT decision function
def sprt_test(residual):
    """ Sequential Probability Ratio Test for anomaly detection """
    likelihood_ratio = stats.norm(mu_1, sigma_0).pdf(residual) / stats.norm(mu_0, sigma_0).pdf(residual)
    return likelihood_ratio > (1 - beta) / alpha
# Apply SPRT
anomalies = np.array([sprt_test(res) for res in residuals])
# Plot results
plt.figure(figsize=(12, 5))
plt.plot(residuals, label="Residuals", color="blue")
plt.axhline(mu_1, color="red", linestyle="dashed", label="Anomaly Threshold")
plt.scatter(np.where(anomalies)[0], residuals[anomalies], color="red", label="Detected Anomalies", zorder=2)
plt.xlabel("Time")
plt.ylabel("Residual Magnitude")
plt.legend()
plt.title("MSET-SPRT Anomaly Detection")
 plt.show()