Descriptive  → What happened?
  "Temperature averaged 27°C last month"

Diagnostic   → Why did it happen?
  "Sensor 3 spiked at 2 AM — the HVAC compressor failed"

Predictive   → What will happen?
  "Based on vibration patterns, motor bearing will fail in ~2 weeks"

Prescriptive → What should we do?
  "Schedule maintenance for motor #4 on Tuesday morning"

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from datetime import datetime, timedelta

# Load sensor data from CSV or InfluxDB export
df = pd.read_csv("sensor_data.csv", parse_dates=["timestamp"])
df = df.set_index("timestamp").sort_index()

# Basic exploration
print(df.shape)                         # (rows, columns)
print(df.dtypes)
print(df.describe())                    # Min, max, mean, std, percentiles
print(df.isnull().sum())               # Missing values per column
print(df.index.min(), df.index.max())  # Date range

# Sample output:
#        temperature  humidity  battery
# count   43200.00   43200.00  43200.00
# mean       24.52      61.20     82.30
# std         2.18       5.80      5.10
# min        18.50      40.00     60.00
# max        38.90      98.00     92.00

# Resample to different time windows
hourly = df.resample('1H').agg({
    'temperature': ['mean', 'min', 'max', 'std'],
    'humidity':    'mean'
})
daily = df.resample('1D').mean()

print(hourly.head())

def zscore_anomalies(series, threshold=3.0, window=None):
    """
    Detect values that are more than `threshold` standard deviations
    from the mean (global or rolling).

    Z-score: z = (x - µ) / σ
    """
    if window:
        # Rolling Z-score (better for non-stationary IoT data)
        rolling_mean = series.rolling(window=window, center=True).mean()
        rolling_std  = series.rolling(window=window, center=True).std()
        z_scores = (series - rolling_mean) / rolling_std
    else:
        # Global Z-score
        z_scores = (series - series.mean()) / series.std()

    anomalies = series[z_scores.abs() > threshold]
    return anomalies, z_scores

df['temperature'] = pd.to_numeric(df['temperature'])

anomalies, z_scores = zscore_anomalies(
    df['temperature'],
    threshold=3.0,
    window=60  # Rolling 60-period window
)

print(f"Found {len(anomalies)} anomalies out of {len(df)} readings")
print(f"Anomaly rate: {len(anomalies)/len(df)*100:.2f}%")
print(anomalies.head(10))

def iqr_anomalies(series, multiplier=1.5):
    """
    Detect outliers using the IQR method (Tukey's fence).
    More robust than Z-score for skewed distributions.

    Lower fence: Q1 - 1.5 * IQR
    Upper fence: Q3 + 1.5 * IQR
    """
    Q1 = series.quantile(0.25)
    Q3 = series.quantile(0.75)
    IQR = Q3 - Q1

    lower_fence = Q1 - multiplier * IQR
    upper_fence = Q3 + multiplier * IQR

    outliers = series[(series < lower_fence) | (series > upper_fence)]

    return outliers, lower_fence, upper_fence

outliers, low, high = iqr_anomalies(df['temperature'], multiplier=1.5)
print(f"Temperature bounds: {low:.1f}°C – {high:.1f}°C")
print(f"Outliers: {len(outliers)}")

def rolling_anomaly_detector(df, column, window=20, threshold=2.5):
    """
    Real-time anomaly detection using a rolling baseline.
    Each new reading is compared against recent history.
    """
    series = df[column].copy()

    roll_mean = series.rolling(window=window).mean()
    roll_std  = series.rolling(window=window).std()

    # Compute confidence bounds
    upper_band = roll_mean + threshold * roll_std
    lower_band = roll_mean - threshold * roll_std

    anomalies = series[(series > upper_band) | (series < lower_band)]

    return {
        'anomalies':   anomalies,
        'upper_band':  upper_band,
        'lower_band':  lower_band,
        'rolling_mean': roll_mean
    }

result = rolling_anomaly_detector(df, 'temperature', window=30, threshold=2.5)
print(f"Rolling anomalies: {len(result['anomalies'])}")

import numpy as np

def analyze_vibration(signal, sample_rate_hz=1000):
    """
    Fast Fourier Transform of vibration signal to detect bearing wear.
    Damaged bearings produce characteristic frequency spikes.
    """
    n = len(signal)
    fft_vals = np.fft.rfft(signal)
    fft_freq  = np.fft.rfffreq(n, d=1.0/sample_rate_hz)

    # Magnitude spectrum
    magnitude = np.abs(fft_vals) / n

    # Find dominant frequencies
    top_n = 5
    peak_indices = np.argsort(magnitude)[-top_n:][::-1]
    peaks = [(fft_freq[i], magnitude[i]) for i in peak_indices]

    print("Top vibration frequencies:")
    for freq, mag in peaks:
        print(f"  {freq:.1f} Hz: magnitude {mag:.4f}")

    # Total harmonic distortion (higher = more wear)
    fundamental_mag = magnitude[np.argmax(magnitude)]
    harmonics_mag   = sum(magnitude[1:]) - fundamental_mag
    thd = harmonics_mag / fundamental_mag if fundamental_mag > 0 else 0

    return {
        "dominant_freq_hz": float(fft_freq[np.argmax(magnitude)]),
        "thd": round(thd, 4),
        "wear_indicator": "HIGH" if thd > 0.5 else "NORMAL"
    }

# Simulate: healthy bearing vs worn bearing
t = np.linspace(0, 1, 1000)
healthy = np.sin(2 * np.pi * 60 * t) + np.random.normal(0, 0.05, 1000)
worn    = np.sin(2 * np.pi * 60 * t) + 0.3 * np.sin(2 * np.pi * 180 * t) + \
          0.2 * np.sin(2 * np.pi * 300 * t) + np.random.normal(0, 0.05, 1000)

print("Healthy bearing:")
print(analyze_vibration(healthy))
print("\nWorn bearing:")
print(analyze_vibration(worn))

def estimate_rul(health_score_history, threshold=0.3):
    """
    Estimate remaining useful life from degrading health scores.
    Uses linear regression on the degradation trend.
    """
    from sklearn.linear_model import LinearRegression
    import numpy as np

    x = np.arange(len(health_score_history)).reshape(-1, 1)
    y = np.array(health_score_history)

    model = LinearRegression()
    model.fit(x, y)

    slope     = model.coef_[0]
    intercept = model.intercept_

    if slope >= 0:  # Improving or stable
        return None, "Stable or improving — no maintenance predicted"

    # Steps until health score hits threshold
    steps_to_threshold = (threshold - intercept) / slope

    return max(0, steps_to_threshold), (
        f"Estimated RUL: {steps_to_threshold:.0f} measurement steps"
    )

# Simulate health degradation
import random
health_scores = [1.0 - i * 0.01 + random.uniform(-0.02, 0.02)
                 for i in range(80)]

rul, message = estimate_rul(health_scores, threshold=0.3)
print(message)
if rul:
    print(f"Recommend maintenance in ~{rul:.0f} cycles")

import paho.mqtt.client as mqtt
import json
import threading
from collections import deque

class StreamAnalyzer:
    """Real-time statistical analysis of incoming MQTT sensor stream."""

    def __init__(self, window=60):
        self._data = {}       # device_id → deque of readings
        self._lock = threading.Lock()
        self._window = window

    def add(self, device_id, value):
        with self._lock:
            if device_id not in self._data:
                self._data[device_id] = deque(maxlen=self._window)
            self._data[device_id].append(value)
            return self._analyze(device_id)

    def _analyze(self, device_id):
        data = list(self._data[device_id])
        if len(data) < 10:
            return None

        mean = sum(data) / len(data)
        std  = (sum((x - mean) ** 2 for x in data) / len(data)) ** 0.5

        latest = data[-1]
        z_score = abs(latest - mean) / std if std > 0 else 0

        return {
            "device_id": device_id,
            "latest":    round(latest, 2),
            "mean":      round(mean, 2),
            "std":       round(std, 2),
            "z_score":   round(z_score, 2),
            "anomaly":   z_score > 3.0,
            "trend":     "UP" if data[-1] > data[-10] else "DOWN"
        }

analyzer = StreamAnalyzer(window=60)

def on_message(client, userdata, msg):
    data = json.loads(msg.payload)
    result = analyzer.add(data['device_id'], data['temperature'])
    if result:
        if result['anomaly']:
            print(f"ANOMALY [{result['device_id']}]: T={result['latest']}°C "
                  f"z={result['z_score']:.1f}")
        else:
            print(f"OK [{result['device_id']}]: T={result['latest']}°C "
                  f"mean={result['mean']}°C trend={result['trend']}")

client = mqtt.Client("stream-analyzer")
client.on_connect = lambda c, u, f, r: c.subscribe("sensors/#", qos=1)
client.on_message = on_message
client.connect("localhost", 1883, 60)
client.loop_forever()