Build a Test Station Health Monitor

Introduction

30 Second Summary

Every factory has test stations that quietly go wrong. Measurements drift, results get flaky, and products slowly creep toward failure thresholds without anyone noticing until customers start complaining.

In this project, you will build a Test Station Health Monitor in a Jupyter notebook that generates realistic production test data, injects three hidden failure modes, and then catches every one of them using statistical and machine learning techniques.

What You'll Build

A Jupyter notebook that acts as a health monitoring dashboard for a simulated production test floor with 6 test stations, using SPC control charts, Kolmogorov-Smirnov drift detection, and Isolation Forest anomaly scoring.

By the end of this project, you'll have:

• A synthetic data generator that models NI TestStand-style production output with three hidden failure modes: drift, flakiness, and silent yield erosion.
• A layered detection pipeline combining SPC control charts, KS drift tests, and Isolation Forest anomaly scoring that catches what traditional pass/fail monitoring misses.
• A unified Station Health Report Card with interactive Plotly visualizations summarizing the health of every test station on your simulated floor.
• A Flaky Station Detector (Secret Mission) that uses consistency scoring and variance analysis to catch intermittently failing stations that game simple yield metrics.

Are there any prerequisites?

You should be comfortable with basic Python and basic pandas (creating and working with DataFrames). You do not need prior experience with statistical process control or machine learning.

Set Up Your Analytics Environment

Before you can detect hidden failures in production test data, you need the right analytical toolkit. This project uses four libraries that each handle a different layer of detection: scipy for statistical drift tests, scikit-learn for machine learning anomaly detection, pycontrolcharts for SPC control charts, and Plotly for interactive visualizations.

In This Step, Get Ready To:

• Install the four Python libraries needed for the project.
• Create a Jupyter notebook and verify all imports load successfully.
• Organize the notebook with section headers for a recruiter-friendly structure.

Install the Required Libraries

Each library in this project serves a specific detection purpose. scipy provides the Kolmogorov-Smirnov test for drift detection. scikit-learn provides Isolation Forest for multivariate anomaly scoring. pycontrolcharts generates SPC control charts with Nelson rule detection. plotly powers the interactive visualizations throughout the notebook.

• Open your terminal and run the following install command.

pip install scikit-learn==1.8.0 scipy==1.17.1 pycontrolcharts[plotly]==0.1.2 plotly==6.7.0

What does this command do?

This installs four libraries with pinned versions so your results are reproducible. The pycontrolcharts[plotly] syntax installs the core library plus its optional Plotly plotting support. You already have pandas and numpy installed as prerequisites, and they also come bundled as dependencies of these libraries.

Why pycontrolcharts instead of building charts manually?

You could calculate control limits and plot charts from scratch using pandas and matplotlib. However, pycontrolcharts handles the heavy lifting of computing XmR statistics and detecting Nelson rule violations in a validated, tested implementation. This lets you focus your time on the anomaly detection analysis instead of chart construction.

Seeing Installation Errors?

pycontrolcharts requires Python 3.10 or higher. Check your Python version by running python --version in your terminal. If you are on an older version, you will need to upgrade before continuing. Help me troubleshoot my installation error.

Create Your Notebook and Verify Imports

Now that the libraries are installed, you need to create the Jupyter notebook where all your analysis will live and confirm everything loads correctly.

• Create a new Jupyter notebook called station_health_monitor.ipynb.
• In the first code cell, paste the following imports and run the cell.

import pandas as pd
import numpy as np
from scipy.stats import ks_2samp
from sklearn.ensemble import IsolationForest
from pycontrolcharts import calc_xmr, plot_control_chart, ChartType
import plotly.express as px
import plotly.graph_objects as go

print("All imports loaded successfully!")

What does each import do?

• pandas and numpy handle data manipulation and numerical operations.
• ks_2samp is the two-sample Kolmogorov-Smirnov test you will use for drift detection.
• IsolationForest is the anomaly detection model that catches multivariate outliers.
• calc_xmr, plot_control_chart, and ChartType are the SPC functions for computing and visualizing control charts.
• plotly.express and plotly.graph_objects power the interactive charts and heatmaps you will build throughout the project.

You should see the message All imports loaded successfully! printed below the cell with no error messages.

I see the success message

All your libraries are installed and loading correctly. You are ready to move on.

I see an error

If you see a ModuleNotFoundError, the library did not install into the same Python environment your Jupyter notebook is using. Try running the install command directly from a notebook cell by adding a ! prefix.

!pip install scikit-learn==1.8.0 scipy==1.17.1 pycontrolcharts[plotly]==0.1.2 plotly==6.7.0

After the install completes, restart your notebook kernel and run the import cell again. Help me fix my import error.

Set Up the Notebook Structure

A well-organized notebook tells a story. Adding section headers upfront gives your analysis a clear narrative arc that is easy for anyone reviewing your work to follow. This is especially important for portfolio projects where a hiring manager might skim the structure before reading the code.

• Add a new Markdown cell at the top of your notebook (above the import cell) with the project title.

# Test Station Health Monitor
Detecting drift, flakiness, and silent yield loss in NI production test data.

• Below your import cell, add five more Markdown cells with the following section headers.

## 1. Data Generation

## 2. SPC Analysis

## 3. Drift Detection

## 4. Anomaly Scoring

## 5. Station Health Report

Why this structure?

Each section maps to a detection layer in your analysis pipeline. You will generate data first, then apply three increasingly sophisticated detection methods (SPC, KS drift tests, Isolation Forest), and finally combine everything into a unified health report. This layered approach mirrors how real test engineering teams build monitoring systems.

Each section header is a placeholder for now. You will add code cells under each section as you work through the remaining steps.

Your analytics environment is set up and your notebook has a clear structure. Next up, you will build a synthetic data generator that creates realistic NI production test data with three hidden failure modes baked in.

Generate Realistic NI Test Station Data

Your health monitor needs something to monitor. In a real NI production environment, TestStand generates thousands of test records per day across multiple stations. Each record captures parametric measurements, pass/fail status, and metadata like station IDs and serial numbers.

Instead of requiring a live SystemLink server (which needs an enterprise license), you will build a synthetic data generator that mirrors the TestStand data model. The real power of this approach is that you can inject controlled failure modes that make the analysis meaningful. Three hidden problems will lurk in your data, and the rest of this project is about catching them.

In This Step, Get Ready To:

• Build a baseline data generator that creates healthy parametric test data for 6 stations over 30 days.
• Inject three hidden failure modes: measurement drift, intermittent flakiness, and silent yield erosion.
• Explore the generated data to confirm it looks normal on the surface.

Build the Baseline Data Generator

NI TestStand structures test data in a hierarchy. A Test Result is the top-level record for one DUT test execution. It contains the program name, status (PASSED/FAILED), serial number, operator, station ID, and timestamps. Inside each result are Test Steps, which hold individual measurements with values, limits, and units.

Your synthetic data flattens this hierarchy into a single DataFrame where each row represents one DUT test result with embedded measurements. This mirrors how production data is typically exported for analysis.

• Add a new code cell in your notebook under the Data Generation section.
• Paste the following function that generates healthy baseline data for 6 stations over 30 days.

import pandas as pd
import numpy as np
from datetime import datetime, timedelta

def generate_baseline_data(n_stations=6, n_days=30, duts_per_station_per_day=100, seed=42):
    """Generate healthy test station data mimicking NI TestStand output."""
    np.random.seed(seed)
    records = []

    # Station-specific baselines (slight natural variation between stations)
    station_baselines = {
        'STATION_01': {'voltage': 5.00, 'current': 2.00, 'resistance': 100.0},
        'STATION_02': {'voltage': 5.02, 'current': 1.98, 'resistance': 100.5},
        'STATION_03': {'voltage': 4.98, 'current': 2.01, 'resistance': 99.8},
        'STATION_04': {'voltage': 5.01, 'current': 1.99, 'resistance': 100.2},
        'STATION_05': {'voltage': 4.99, 'current': 2.02, 'resistance': 100.1},
        'STATION_06': {'voltage': 5.00, 'current': 2.00, 'resistance': 99.9},
    }

    # Spec limits for pass/fail determination
    spec_limits = {
        'voltage':    (4.5, 5.5),
        'current':    (1.5, 2.5),
        'resistance': (90.0, 110.0),
    }

    # Measurement noise (standard deviation)
    noise_std = {'voltage': 0.05, 'current': 0.03, 'resistance': 1.0}

    start_date = datetime(2025, 6, 1)
    serial_counter = 0

    for day in range(n_days):
        timestamp_base = start_date + timedelta(days=day)
        for station_id, baselines in station_baselines.items():
            for dut in range(duts_per_station_per_day):
                serial_counter += 1
                # Spread DUTs across the workday (8 hours)
                timestamp = timestamp_base + timedelta(
                    hours=8 + (dut / duts_per_station_per_day) * 8,
                    minutes=np.random.randint(0, 60)
                )

                # Generate measurements from normal distributions
                voltage = np.random.normal(baselines['voltage'], noise_std['voltage'])
                current = np.random.normal(baselines['current'], noise_std['current'])
                resistance = np.random.normal(baselines['resistance'], noise_std['resistance'])

                # Determine pass/fail based on spec limits
                passed = (
                    spec_limits['voltage'][0] <= voltage <= spec_limits['voltage'][1]
                    and spec_limits['current'][0] <= current <= spec_limits['current'][1]
                    and spec_limits['resistance'][0] <= resistance <= spec_limits['resistance'][1]
                )

                records.append({
                    'timestamp': timestamp,
                    'station_id': station_id,
                    'serial_number': f'DUT-{serial_counter:06d}',
                    'voltage': round(voltage, 4),
                    'current': round(current, 4),
                    'resistance': round(resistance, 4),
                    'status': 'PASSED' if passed else 'FAILED',
                    'day': day + 1,
                })

    return pd.DataFrame(records)

# Generate the baseline (healthy) dataset
df = generate_baseline_data()
print(f"Generated {len(df):,} records")
print(f"Stations: {df['station_id'].nunique()}")
print(f"Date range: {df['timestamp'].min().date()} to {df['timestamp'].max().date()}")

What does this code do?

This function creates 18,000 test records (6 stations x 30 days x 100 DUTs per station per day). Each record includes three parametric measurements drawn from normal distributions with station-specific baselines.

The spec limits define the pass/fail boundaries. Voltage must be between 4.5V and 5.5V, current between 1.5A and 2.5A, and resistance between 90 and 110 ohms. A DUT passes only if all three measurements fall within their respective limits.

Why synthetic data instead of real NI SystemLink data?

NI SystemLink is enterprise software that requires a license and a running server. Synthetic data lets you work for free while injecting controlled failure modes that make the analysis teachable. If you have access to real NI data, you can swap it into this same pipeline later.

Inject Three Hidden Failure Modes

Healthy data is not very interesting to analyze. Real production test floors develop subtle problems that traditional pass/fail monitoring completely misses. You will now inject three failure patterns into the data, each representing a different category of silent failure.

• Add a new code cell below your baseline generator.
• Paste the following function that modifies the clean data to inject the three failure modes.

def inject_failure_modes(df):
    """Inject three hidden failure modes into the dataset."""
    df = df.copy()

    spec_limits = {
        'voltage':    (4.5, 5.5),
        'current':    (1.5, 2.5),
        'resistance': (90.0, 110.0),
    }

    # --- Failure Mode 1: Gradual drift on Station 3 ---
    # Voltage mean shifts from ~5.0V toward 5.35V over days 15-30
    mask_s3 = (df['station_id'] == 'STATION_03') & (df['day'] >= 15)
    drift_days = df.loc[mask_s3, 'day'] - 15  # 0 to 15
    drift_amount = (drift_days / 15) * 0.35   # Linear drift up to +0.35V
    df.loc[mask_s3, 'voltage'] += drift_amount.values

    # --- Failure Mode 2: Intermittent flakiness on Station 5 ---
    # Random 10% of measurements get noise spikes (3-5x normal std dev)
    mask_s5 = df['station_id'] == 'STATION_05'
    n_s5 = mask_s5.sum()
    flaky_mask = np.random.random(n_s5) < 0.10  # 10% of Station 5 DUTs
    flaky_indices = df.loc[mask_s5].index[flaky_mask]
    for param, std in [('voltage', 0.05), ('current', 0.03), ('resistance', 1.0)]:
        spike_multiplier = np.random.uniform(3, 5, size=len(flaky_indices))
        spike_direction = np.random.choice([-1, 1], size=len(flaky_indices))
        df.loc[flaky_indices, param] += spike_direction * spike_multiplier * std

    # --- Failure Mode 3: Silent yield erosion on Station 2 ---
    # All parameters slowly creep toward upper spec limits over 30 days
    mask_s2 = df['station_id'] == 'STATION_02'
    creep_progress = (df.loc[mask_s2, 'day'] - 1) / 29  # 0.0 to 1.0
    df.loc[mask_s2, 'voltage'] += (creep_progress * 0.28).values
    df.loc[mask_s2, 'current'] += (creep_progress * 0.30).values
    df.loc[mask_s2, 'resistance'] += (creep_progress * 6.5).values

    # Recalculate pass/fail status after injecting failures
    df['status'] = df.apply(
        lambda row: 'PASSED' if (
            spec_limits['voltage'][0] <= row['voltage'] <= spec_limits['voltage'][1]
            and spec_limits['current'][0] <= row['current'] <= spec_limits['current'][1]
            and spec_limits['resistance'][0] <= row['resistance'] <= spec_limits['resistance'][1]
        ) else 'FAILED',
        axis=1
    )

    return df

# Inject failure modes into the baseline data
df = inject_failure_modes(df)
print("Failure modes injected.")
print(f"Overall yield: {(df['status'] == 'PASSED').mean():.1%}")
print(f"\nYield by station:")
print(df.groupby('station_id')['status'].apply(lambda x: (x == 'PASSED').mean()).to_string())

What does this code do?

This function modifies the clean dataset in three targeted ways, then recalculates pass/fail status so the labels stay consistent with the new measurements.

• Station 3 drift: voltage mean shifts from ~5.0V toward 5.35V linearly over days 15-30. It stays within the 5.5V spec limit but trends dangerously close.
• Station 5 flakiness: 10% of DUTs randomly receive noise spikes 3-5x the normal standard deviation across all parameters. Some will trip pass/fail, but overall yield stays above 90%.
• Station 2 silent yield erosion: all three parameters slowly creep toward their upper spec limits over the full 30 days. No single measurement triggers a failure, but the correlated shift is suspicious.

Why are these failure modes dangerous in real production?

Drift causes out-of-spec escapes. A station that slowly drifts will eventually ship bad parts, but traditional yield dashboards will not flag it until failures spike.

Flakiness masks real defects. If a station intermittently fails DUTs that are actually good, operators learn to re-test and ignore failures. This trains the team to dismiss real problems.

Silent yield erosion leads to field failures. Parts that technically pass but sit near spec limits are fragile. They fail under real-world conditions that the test floor does not simulate.

Explore the Generated Data

Before building any detection logic, you should look at the data the way a test engineer would on a Monday morning. Summary statistics and a quick time series plot. The key insight here is that everything looks fine at first glance.

• Add a new code cell and paste the following exploration code.

# Check the shape and structure
print(f"Dataset shape: {df.shape}")
print(f"\nRecords per station:")
print(df['station_id'].value_counts().sort_index().to_string())

# Summary statistics grouped by station
print(f"\nVoltage stats by station:")
print(df.groupby('station_id')['voltage'].describe().round(4).to_string())

print(f"\nCurrent stats by station:")
print(df.groupby('station_id')['current'].describe().round(4).to_string())

print(f"\nResistance stats by station:")
print(df.groupby('station_id')['resistance'].describe().round(4).to_string())

The summary statistics might show slightly elevated means for Station 2 and Station 3, but nothing that screams "problem" in a table of numbers. This is exactly how silent failures hide in real production data.

• Add another code cell to create a time series plot of voltage measurements per station.

import plotly.express as px

# Compute daily average voltage per station
daily_voltage = df.groupby(['day', 'station_id'])['voltage'].mean().reset_index()

fig = px.line(
    daily_voltage,
    x='day',
    y='voltage',
    color='station_id',
    title='Daily Average Voltage by Station',
    labels={'day': 'Day', 'voltage': 'Voltage (V)', 'station_id': 'Station'},
)

# Add spec limit reference lines
fig.add_hline(y=5.5, line_dash='dash', line_color='red', annotation_text='Upper Spec Limit')
fig.add_hline(y=4.5, line_dash='dash', line_color='red', annotation_text='Lower Spec Limit')

fig.update_layout(height=500, template='plotly_white')
fig.show()

What should I see?

The plot shows daily average voltage for each station. Most stations hold steady around 5.0V. Station 3 starts climbing around day 15, trending upward toward the upper spec limit. Station 2 shows a gentler upward creep. All lines stay within the red spec limit boundaries.

This is the problem your health monitor will solve. A quick glance at this chart might catch Station 3 if you are paying close attention. Station 2's creep is much harder to spot. Station 5's flakiness does not show up at all in daily averages because the spikes average out.

Your dataset is ready with 18,000 records and three hidden problems waiting to be caught. Next up, you will apply Statistical Process Control charts to start peeling back the layers.

Build SPC Control Charts per Station

Your synthetic dataset is ready, and a quick glance at the raw data suggests everything is running smoothly. That is exactly the problem. To catch hidden process issues, you need a technique that separates normal variation from signals of trouble. Statistical Process Control (SPC) using XmR charts is the foundational technique that test engineers use to monitor individual measurements against natural process limits. This is your first detection layer.

In This Step, Get Ready To:

• Calculate XmR control charts for each station's voltage measurements using pycontrolcharts.
• Visualize control charts and identify Nelson rule violations that flag process instability.
• Summarize SPC findings across all stations and identify which failure modes SPC catches and which it misses.

Calculate XmR Charts for Each Station

An XmR chart (Individuals and Moving Range) plots each individual measurement against control limits calculated from the data's own natural variation. Unlike X-bar charts that require subgroups of measurements, XmR charts work on individual readings. This makes them the right choice here because each DUT produces one measurement per parameter, not a batch.

The pycontrolcharts library handles the heavy lifting. Its calc_xmr() function takes a list of numeric values and returns a pandas DataFrame with the center line, upper and lower control limits, and any Nelson rule violations already flagged.

• Add a new cell in your notebook under the SPC Analysis section.
• Paste the following code to compute XmR chart data for the voltage parameter across all 6 stations.

from pycontrolcharts import calc_xmr

# Compute XmR chart data for voltage measurements, one chart per station
xmr_results = {}

for station in sorted(df['station_id'].unique()):
    station_data = df[df['station_id'] == station].sort_values('timestamp')
    voltage_values = station_data['voltage'].tolist()

    # calc_xmr returns a DataFrame with control limits and violation flags
    xmr_df = calc_xmr(voltage_values)
    xmr_results[station] = xmr_df

    # Count violations for a quick summary
    violation_count = xmr_df[xmr_df['violations'].apply(len) > 0].shape[0]
    print(f"{station}: {violation_count} Nelson rule violations detected")

What does this code do?

For each station, you extract the voltage measurements in chronological order and pass them to calc_xmr(). The returned DataFrame includes columns like value, center_line, ucl (upper control limit), lcl (lower control limit), and violations. The violations column contains a list of Nelson rule codes for each data point. An empty list means the point is in control.

You should see Station 3 (the drifting station) report a noticeably higher violation count than the healthy stations. Stations 1, 4, and 6 should have few or zero violations.

Visualize Control Charts and Identify Violations

Numbers alone do not tell the full story. Plotting the control charts side by side for a healthy station and a drifting station makes the difference visually obvious.

• Add a new cell and paste the following code to plot control charts for Station 1 (healthy) and Station 3 (drifting).

from pycontrolcharts import plot_control_chart, ChartType

# Plot a healthy station's control chart
fig_healthy = plot_control_chart(
    xmr_results['STATION_01'],
    ChartType.XMR,
    title="Station 01 - Voltage (Healthy)"
)
fig_healthy.show()

# Plot the drifting station's control chart
fig_drift = plot_control_chart(
    xmr_results['STATION_03'],
    ChartType.XMR,
    title="Station 03 - Voltage (Drifting)"
)
fig_drift.show()

What should I see?

Station 01's chart should show measurements scattered randomly around the center line, staying well within the upper and lower control limits. Station 03's chart should show a clear upward trend starting around the midpoint, with data points clustering near or beyond the upper control limit in the later half.

Now extract the specific Nelson rule violations from Station 3's chart data to understand what patterns the SPC analysis detected.

• Add a new cell and paste the following code to filter and display the violations.

# Extract rows with Nelson rule violations for Station 03
station3_xmr = xmr_results['STATION_03']
violations_df = station3_xmr[station3_xmr['violations'].apply(len) > 0]

print(f"Station 03 has {len(violations_df)} points with violations:\n")
for idx, row in violations_df.iterrows():
    print(f"  Point {row['point_id']}: value={row['value']:.4f}, violations={row['violations']}")

What Are Nelson Rules?

Nelson rules are a set of 8 decision rules that detect non-random patterns in control chart data. The most relevant ones for this project are:

• Rule 1: A single point beyond 3-sigma (an extreme outlier).
• Rule 2: Nine consecutive points on the same side of the center line (a sustained shift).
• Rule 3: Six consecutive points steadily increasing or decreasing (a trend).
• Rule 6: Four of five consecutive points beyond 1-sigma on the same side (a small but persistent shift).

Station 3's drift should trigger Rules 2, 3, and 6 because the voltage measurements are trending upward and clustering on one side of the center line.

Summarize SPC Findings Across All Stations

You have looked at individual charts, but a test floor has multiple stations and multiple parameters. You need a summary view that ranks stations by their SPC health so you can quickly spot which ones need attention.

• Add a new cell and paste the following code to build a summary DataFrame of violations across all stations and all three parameters.

# Build SPC summary across all stations and all parameters
parameters = ['voltage', 'current', 'resistance']
spc_summary = []

for station in sorted(df['station_id'].unique()):
    station_data = df[df['station_id'] == station].sort_values('timestamp')
    row = {'station_id': station}

    for param in parameters:
        values = station_data[param].tolist()
        xmr_df = calc_xmr(values)
        violation_count = xmr_df[xmr_df['violations'].apply(len) > 0].shape[0]
        row[f'{param}_violations'] = violation_count

    # Total violations across all parameters
    row['total_violations'] = sum(row[f'{p}_violations'] for p in parameters)
    spc_summary.append(row)

spc_summary_df = pd.DataFrame(spc_summary).sort_values('total_violations', ascending=False)
print(spc_summary_df.to_string(index=False))

What does this code do?

This loops through every station-parameter combination (6 stations x 3 parameters = 18 XmR charts), counts the Nelson rule violations for each, and assembles a summary table. The total_violations column gives you a quick ranking of which stations are most unstable.

Station 3 should rank highest in total violations due to its voltage drift. The critical insight here is what SPC does not catch. Station 2 (silent yield erosion) may show only a few violations because its parameters are creeping slowly and uniformly. Station 5 (flaky) may show scattered violations but not a clear pattern. This is exactly why you need additional detection layers.

Pro tip: save your SPC results

Store spc_summary_df as a variable you can reference later. You will combine it with drift detection and anomaly scoring results in Step 6 to build the unified Station Health Report Card.

Your SPC analysis has surfaced Station 3's drift loud and clear. But the summary also revealed a gap: slow, correlated parameter creep and intermittent noise do not always trigger Nelson rules. Next up, you will use the Kolmogorov-Smirnov test to compare measurement distributions against a golden baseline and catch the drift that SPC missed.

Detect Measurement Drift with KS Tests

Your SPC control charts flagged Station 3's obvious drift, but what about Station 2's slow, silent creep toward spec limits? Control charts excel at catching individual outlier points and runs, but they can miss gradual distributional shifts where the entire measurement population slides sideways.

The Kolmogorov-Smirnov (KS) test solves this by comparing two entire distributions and asking: "Did the shape of this station's measurements change?" In this step, you'll define a golden baseline from when all stations were healthy, then compare recent measurements against it to quantify exactly how much each station has drifted.

In This Step, Get Ready To:

• Define a golden baseline window and a recent comparison window from your test data.
• Run two-sample KS tests for every station-parameter combination to quantify drift.
• Build a drift severity heatmap that reveals which stations are shifting away from their known-good state.

Define a Golden Baseline and Comparison Windows

In NI production testing, a golden baseline represents the known-good operating state of a station right after calibration. Every measurement taken during this window is trusted. By comparing recent measurements against this baseline, you can detect whether a station's behavior has changed, even if every individual measurement still passes spec limits.

You need two time slices from your dataset: a baseline window from when all stations were healthy, and a recent window from after the failure modes have had time to develop.

• Add a new markdown cell in your notebook titled Drift Detection.
• In a new code cell, slice your DataFrame into the two windows.

# Golden baseline: days 1-7 (all stations healthy after calibration)
baseline = df[df['day'] <= 7]

# Recent window: days 23-30 (failure modes fully developed)
recent = df[df['day'] >= 23]

print(f"Baseline records: {len(baseline)}")
print(f"Recent records: {len(recent)}")
print(f"Stations in baseline: {baseline['station_id'].nunique()}")
print(f"Stations in recent: {recent['station_id'].nunique()}")

Why Days 1-7 and 23-30?

Days 1-7 are before any injected failure modes begin (Station 3's drift starts at day 15, and Station 2's erosion is minimal in the first week). Days 23-30 give the failure modes enough time to develop measurable shifts. This gap between windows makes the comparison more meaningful.

Run KS Tests per Station per Parameter

The two-sample KS test measures the maximum difference between two cumulative distribution functions (CDFs). A large KS statistic means the two distributions are shaped differently. The test also returns a p-value. A small p-value (below 0.05) means the distributions are statistically different, which signals drift.

You'll run this test for every combination of station and parameter (voltage, current, resistance), producing 18 comparisons total.

• Add the following code to loop through all station-parameter pairs and collect KS results.

from scipy.stats import ks_2samp

stations = df['station_id'].unique()
parameters = ['voltage', 'current', 'resistance']

ks_results = []

for station in stations:
    for param in parameters:
        # Get baseline and recent measurements for this station-parameter pair
        baseline_values = baseline[baseline['station_id'] == station][param]
        recent_values = recent[recent['station_id'] == station][param]

        # Run the two-sample KS test
        stat, p_value = ks_2samp(baseline_values, recent_values)

        ks_results.append({
            'station_id': station,
            'parameter': param,
            'ks_statistic': round(stat, 4),
            'p_value': p_value,
            'drift_detected': p_value < 0.05
        })

ks_df = pd.DataFrame(ks_results)
print(ks_df.to_string(index=False))

What Does This Code Do?

For each of the 6 stations and 3 parameters, the code extracts the baseline measurements (days 1-7) and recent measurements (days 23-30). It passes both arrays to ks_2samp(), which returns the KS statistic (maximum CDF difference) and a p-value.

The drift_detected column flags any pair where the p-value is below 0.05, meaning there is statistically significant evidence that the distribution has changed.

What Should I See?

Station 3 should show drift_detected: True for voltage with a very small p-value. Station 2 should show moderate drift across multiple parameters (voltage, current, and resistance all creeping). Healthy stations like Station 1, 4, and 6 should show drift_detected: False with p-values well above 0.05.

Build a Drift Severity Heatmap

A table of p-values is useful, but a heatmap makes the drift patterns immediately visible. You'll plot -log10(p-value) for each station-parameter pair. This transformation makes the visualization intuitive: higher values mean stronger evidence of drift. A value above ~1.3 corresponds to p < 0.05 (statistically significant).

• Add the following code to create the heatmap.

import plotly.express as px
import numpy as np

# Calculate -log10(p-value) for visualization
# Clip p-values to avoid log(0)
ks_df['neg_log_p'] = -np.log10(ks_df['p_value'].clip(lower=1e-300))

# Pivot into a matrix for the heatmap
heatmap_data = ks_df.pivot(
    index='station_id',
    columns='parameter',
    values='neg_log_p'
)

fig = px.imshow(
    heatmap_data,
    labels=dict(x='Parameter', y='Station', color='-log10(p-value)'),
    title='Drift Severity Heatmap: Station vs Parameter',
    color_continuous_scale='RdYlGn_r',
    aspect='auto',
    text_auto='.1f'
)

# Add a reference line annotation for significance threshold
fig.update_layout(
    coloraxis_colorbar_title='-log10(p)',
    width=700,
    height=450
)

fig.show()

What Does This Heatmap Tell You?

The RdYlGn_r color scale maps green to low drift (healthy) and red to high drift (problematic). Station 3 should light up red for voltage, confirming the drift your SPC charts already caught.

The more interesting finding is Station 2. It should show moderate orange/red across all three parameters. This is the silent yield erosion pattern. No single parameter triggers an alarm on its own, but the correlated creep across voltage, current, and resistance is clearly visible in the heatmap. SPC alone would likely miss this because each individual measurement stays within control limits.

Why KS Tests Instead of Other Drift Methods?

Alternatives like Population Stability Index (PSI) require binning your data, which introduces artifacts. Jensen-Shannon divergence is more information-theoretic and less intuitive for manufacturing engineers. The KS test is non-parametric (works on any continuous distribution), returns both a statistic and p-value for clear interpretation, and is the industry standard for drift detection in manufacturing.

Your KS tests confirmed Station 3's drift and uncovered Station 2's silent multi-parameter creep. But both of these methods analyze one parameter at a time. Next, you'll use Isolation Forest to catch multivariate anomalies where the combination of parameter values is suspicious, even when each individual value looks fine.

Score Anomalies with Isolation Forest

Your SPC charts caught Station 3's drift, and your KS tests revealed Station 2's silent creep. But both of those techniques analyze one parameter at a time. What if a DUT's voltage, current, and resistance all shift together in a correlated way that no single-parameter check would flag? That pattern requires multivariate anomaly detection. In this step, you'll use Isolation Forest to score every DUT across all parameters simultaneously, catching suspicious "passing" units that your earlier layers missed.

In This Step, Get Ready To:

• Build a multivariate feature matrix with engineered features for each DUT.
• Train an Isolation Forest on healthy baseline data and score all DUTs for anomalies.
• Visualize anomaly clusters by station and compare results with SPC and KS findings.

Prepare the Multivariate Feature Matrix

Single-parameter analysis treats each measurement independently. A DUT where voltage is slightly high, current is slightly high, and resistance is slightly high might look fine on each individual chart. But the fact that all three shifted together is a strong signal that something is wrong with the station or the DUT. To capture these correlated shifts, you need to analyze all parameters as a single feature vector per DUT.

You'll also engineer additional features that make anomalies easier to detect. A raw voltage of 5.2V means different things depending on the station's typical average. Converting measurements to z-scores and deviations from station means gives the model a normalized view.

Add a new code cell under your Anomaly Scoring section:

# Build feature matrix for Isolation Forest
# Each row = one DUT, columns = raw measurements + engineered features

# Compute per-station means and stds from the full dataset
station_stats = df.groupby('station_id')[['voltage', 'current', 'resistance']].agg(['mean', 'std'])

def engineer_features(row):
    station = row['station_id']
    features = {}

    for param in ['voltage', 'current', 'resistance']:
        # Raw measurement
        features[param] = row[param]

        # Deviation from station mean
        st_mean = station_stats.loc[station, (param, 'mean')]
        st_std = station_stats.loc[station, (param, 'std')]

        features[f'{param}_deviation'] = row[param] - st_mean

        # Z-score per parameter
        features[f'{param}_zscore'] = (row[param] - st_mean) / st_std if st_std > 0 else 0

    return pd.Series(features)

# Apply feature engineering to every DUT record
feature_matrix = df.apply(engineer_features, axis=1)

print(f"Feature matrix shape: {feature_matrix.shape}")
print(f"Columns: {list(feature_matrix.columns)}")
feature_matrix.head()

What does this code do?

For each DUT record, this creates 9 features: the 3 raw measurements (voltage, current, resistance), 3 deviation-from-station-mean values, and 3 z-scores. The deviations capture how far each measurement is from what's typical for that specific station. The z-scores normalize everything to a common scale so that voltage (around 5V) and resistance (around 100 ohms) are comparable.

Why multivariate analysis matters

Station 2's silent yield erosion shifts all three parameters toward their spec limits simultaneously. Each individual shift is small enough to pass single-parameter checks. But when you look at all three together, the correlated movement stands out. This is exactly the kind of pattern that Isolation Forest is designed to catch.

Train Isolation Forest and Score DUTs

Isolation Forest works on a simple but powerful principle. It builds random decision trees that partition data. Anomalies are data points that are easy to isolate, meaning they get separated from the rest of the data in fewer splits. Normal points are buried deep in dense clusters and require many splits to isolate. The model assigns a score based on the average path length across all trees.

You'll train the model only on baseline data from days 1-7, when all stations were healthy. This teaches the model what "normal" looks like. Then you'll score every DUT in the dataset to find which ones deviate from that learned normal.

• Add a new code cell with the following:

from sklearn.ensemble import IsolationForest

# Split feature matrix into baseline (days 1-7) and full dataset
baseline_mask = df['day'] <= 7
baseline_features = feature_matrix[baseline_mask]

print(f"Baseline samples (days 1-7): {len(baseline_features)}")
print(f"Total samples: {len(feature_matrix)}")

# Train Isolation Forest on healthy baseline data
iso_forest = IsolationForest(
    contamination=0.05,   # Expect ~5% of data to be anomalous
    random_state=42,      # Reproducible results
    n_estimators=100      # Number of isolation trees
)
iso_forest.fit(baseline_features)

# Score ALL DUTs (baseline + post-baseline)
df['anomaly_score'] = iso_forest.score_samples(feature_matrix)
df['anomaly_label'] = iso_forest.predict(feature_matrix)

# Count anomalies
n_anomalies = (df['anomaly_label'] == -1).sum()
n_total = len(df)
print(f"\nAnomalies detected: {n_anomalies} / {n_total} ({n_anomalies/n_total*100:.1f}%)")
print(f"\nAnomalies per station:")
print(df[df['anomaly_label'] == -1].groupby('station_id').size())

What is contamination in Isolation Forest?

The contamination parameter tells the model what fraction of the training data you expect to be anomalous. Setting it to 0.05 means roughly 5%. The model uses this to set the decision threshold for labeling points as normal (1) or anomalous (-1). The score_samples() method returns continuous anomaly scores where lower values mean more anomalous.

What should I see?

You should see that the anomaly count is concentrated on Stations 2, 3, and 5. Station 3 (drift) and Station 5 (flakiness) should have noticeable anomaly counts. Station 2 (silent yield erosion) is the key finding here. Its DUTs technically pass, but the correlated parameter shifts make them stand out to the multivariate model.

Visualize Anomaly Clusters by Station

Numbers alone don't tell the full story. A scatter plot of anomaly scores over time, colored by station, reveals the patterns behind the numbers. You'll see how each failure mode creates a distinct visual signature.

• Add a new code cell to create the anomaly visualization:

import plotly.express as px

# Create scatter plot of anomaly scores over time
fig = px.scatter(
    df,
    x='day',
    y='anomaly_score',
    color='station_id',
    title='Isolation Forest Anomaly Scores by Station Over Time',
    labels={
        'day': 'Day',
        'anomaly_score': 'Anomaly Score (lower = more anomalous)',
        'station_id': 'Station'
    },
    opacity=0.5
)

# Add threshold line
threshold = df[df['anomaly_label'] == -1]['anomaly_score'].max()
fig.add_hline(
    y=threshold,
    line_dash='dash',
    line_color='red',
    annotation_text='Anomaly Threshold'
)

fig.update_layout(height=500)
fig.show()

# Show anomaly rate by station and time period
print("\nAnomaly rate by station (days 1-7 vs days 23-30):")
for station in sorted(df['station_id'].unique()):
    early = df[(df['station_id'] == station) & (df['day'] <= 7)]
    late = df[(df['station_id'] == station) & (df['day'] >= 23)]
    early_rate = (early['anomaly_label'] == -1).mean() * 100
    late_rate = (late['anomaly_label'] == -1).mean() * 100
    print(f"  {station}: {early_rate:.1f}% (early) -> {late_rate:.1f}% (late)")

What should I see?

The scatter plot should reveal three distinct patterns. Station 3's anomalies cluster in the later days (15-30) as its drift intensifies. Station 5's anomalies appear as scattered spikes throughout the timeline, reflecting its intermittent flakiness. Station 2's anomalies gradually increase in the later days as all three parameters creep toward spec limits together. Healthy stations (1, 4, 6) should have very few points below the threshold.

The anomaly rate comparison between early and late periods is especially telling. Healthy stations should show similar rates in both windows. Problem stations should show a clear jump in their late-period anomaly rate.

How the Three Detection Layers Complement Each Other

Each layer catches different things. SPC control charts excel at catching Station 3's drift because Nelson rules are designed for trending patterns. KS tests catch Station 2's distributional shift by comparing entire distributions rather than individual points. Isolation Forest catches the multivariate correlation in Station 2's data that single-parameter methods miss entirely.

No single technique catches all three failure modes. That's why production test monitoring systems layer multiple detection methods. In the next step, you'll combine all three signals into a unified health score.

You now have three independent detection layers, each revealing different aspects of station health. Next up, you'll combine all of these signals into a unified Station Health Report Card that gives a single, actionable view of your entire test floor.

Build the Unified Station Health Report

You now have three independent detection layers, each catching different failure patterns. SPC flagged Station 3's drift. KS tests revealed Station 2's silent creep. Isolation Forest surfaced multivariate anomalies across multiple stations. The problem is that no single method caught everything. A test engineer checking only one of these signals would miss critical issues.

In this step, you will combine all three detection signals into a single Station Health Report Card. This gives you one place to look for the overall health of your test floor, with clear, actionable status labels.

In This Step, Get Ready To:

• Compute a composite health score per station from SPC violations, KS drift magnitude, and anomaly rate.
• Generate a Station Health Report Card with status labels and a Plotly bar chart.
• Write a narrative interpretation of the results like a real test engineering report.

Compute a Composite Health Score per Station

Each of your three detection methods produces a different type of output. SPC gives you violation counts. KS tests give you p-values. Isolation Forest gives you anomaly labels. To compare stations fairly, you need to normalize each signal to a common 0-100 scale and then combine them with weights that reflect their importance.

• Add a new markdown cell with the heading ## Station Health Report.
• In a new code cell, create a scoring function that takes in the SPC summary, KS results, and Isolation Forest predictions, then computes a composite health score per station.

def compute_health_scores(spc_summary, ks_results, df_scored):
    """Combine SPC, KS, and Isolation Forest signals into a 0-100 health score."""
    health_data = []

    for station in df_scored['station_id'].unique():
        # SPC signal: total violation count for this station
        station_spc = spc_summary[spc_summary['station_id'] == station]
        spc_violations = station_spc['violation_count'].sum() if len(station_spc) > 0 else 0

        # KS signal: average -log10(p-value) across parameters
        station_ks = ks_results[ks_results['station_id'] == station]
        avg_drift = station_ks['neg_log10_pvalue'].mean() if len(station_ks) > 0 else 0

        # Isolation Forest signal: percentage of DUTs flagged as anomalies
        station_data = df_scored[df_scored['station_id'] == station]
        anomaly_rate = (station_data['anomaly_label'] == -1).mean() * 100

        health_data.append({
            'station_id': station,
            'spc_violations': spc_violations,
            'drift_severity': round(avg_drift, 2),
            'anomaly_rate_pct': round(anomaly_rate, 2)
        })

    health_df = pd.DataFrame(health_data)

    # Normalize each signal to 0-100 (higher = worse)
    for col in ['spc_violations', 'drift_severity', 'anomaly_rate_pct']:
        col_max = health_df[col].max()
        if col_max > 0:
            health_df[f'{col}_norm'] = (health_df[col] / col_max) * 100
        else:
            health_df[f'{col}_norm'] = 0

    # Weighted composite: SPC 30%, KS drift 40%, Isolation Forest 30%
    health_df['risk_score'] = (
        health_df['spc_violations_norm'] * 0.30 +
        health_df['drift_severity_norm'] * 0.40 +
        health_df['anomaly_rate_pct_norm'] * 0.30
    ).round(1)

    # Health score is inverse of risk (100 = perfectly healthy)
    health_df['health_score'] = (100 - health_df['risk_score']).round(1)

    return health_df

health_df = compute_health_scores(spc_summary, ks_results, df_scored)
health_df[['station_id', 'spc_violations', 'drift_severity', 'anomaly_rate_pct', 'health_score']]

Why These Weights?

The weighting gives KS drift the highest influence (40%) because distributional shift is the most dangerous failure mode in production. It causes out-of-spec escapes that reach customers. SPC and Isolation Forest each get 30% because they catch complementary patterns: obvious outliers and subtle multivariate correlations.

In a real deployment, you would tune these weights based on historical data about which failure modes caused the most costly escapes.

Generate the Station Health Report Card

A numeric score is useful for ranking, but a test engineer glancing at a dashboard needs an instant read. You will add status labels that translate scores into actionable categories and then visualize the results.

• In a new code cell, add status labels and display the final report card.

# Assign status labels based on health score thresholds
def assign_status(score):
    if score >= 80:
        return 'Healthy'
    elif score >= 50:
        return 'Warning'
    else:
        return 'Critical'

health_df['status'] = health_df['health_score'].apply(assign_status)

# Build the final report card
report_card = health_df[[
    'station_id', 'health_score', 'spc_violations',
    'drift_severity', 'anomaly_rate_pct', 'status'
]].sort_values('health_score', ascending=True)

print("=== Station Health Report Card ===")
report_card

What Should I See?

Your report card should show Station 3 and Station 2 with the lowest health scores (Warning or Critical status). Station 5 should also show elevated anomaly rates from its flaky measurements. The remaining healthy stations (1, 4, 6) should score above 80 with Healthy status.

• In a new code cell, create a Plotly bar chart to visualize health scores across stations.

import plotly.express as px

# Color-code bars by status
color_map = {'Healthy': '#2ecc71', 'Warning': '#f39c12', 'Critical': '#e74c3c'}

fig = px.bar(
    report_card,
    x='station_id',
    y='health_score',
    color='status',
    color_discrete_map=color_map,
    title='Station Health Report Card',
    labels={'health_score': 'Health Score (0-100)', 'station_id': 'Station'},
    text='health_score'
)

fig.update_layout(
    yaxis_range=[0, 105],
    showlegend=True
)
fig.update_traces(textposition='outside')
fig.show()

What Does This Chart Show?

The bar chart gives you an at-a-glance view of your entire test floor. Green bars mean healthy stations that need no attention. Yellow bars are warnings that should be investigated during the next maintenance window. Red bars are critical stations that need immediate action.

Add Narrative Interpretation

Numbers and charts tell part of the story. A real test engineering report includes written interpretation that explains what the data means and what actions to take. This is what separates a data dump from actionable intelligence.

• Add a new markdown cell below your bar chart with a report-style interpretation of the results.

### Interpretation & Recommended Actions

**Station 3 (Voltage Drift):** Shows clear upward drift in voltage measurements
starting around day 15. SPC control charts flagged multiple Nelson rule violations.
KS test confirms statistically significant distributional shift.
**Action:** Schedule recalibration of voltage measurement fixture.

**Station 2 (Silent Yield Erosion):** All three parameters are slowly creeping
toward their upper spec limits. No individual measurement has tripped a failure,
but the correlated shift across voltage, current, and resistance suggests a
systematic issue. Isolation Forest flagged multiple "passing" DUTs as anomalous.
**Action:** Investigate shared fixture components and environmental factors.

**Station 5 (Intermittent Flakiness):** Elevated anomaly rate from sporadic
measurement spikes. Overall yield appears acceptable, but individual DUTs
experience random failures.
**Action:** Inspect probe contacts and cable connections for intermittent faults.

**Stations 1, 4, 6:** Operating within normal parameters across all detection
methods. No action required.

How Would This Feed Into Real NI Workflows?

In a real NI production environment, this type of analysis would feed into NI SystemLink dashboards. SystemLink aggregates test data from TestStand stations across the factory floor. The health scores and status labels you built here could be published as custom dashboard widgets, giving test engineers a live view of station health.

The recommended actions would become work orders in a maintenance tracking system, closing the loop between detection and resolution.

Pro tip: make it portfolio-ready

Add a final markdown cell at the top of your notebook with a brief project summary. Include what problem it solves, what techniques you used, and what you found. Recruiters and hiring managers often skim the first cell to decide whether to keep reading.

Your Station Health Report Card is complete. You have built a layered detection system that catches drift, silent yield erosion, and flakiness, then rolls everything into a single actionable view. Next up, take on the Secret Mission to build a dedicated Flaky Station Detector that goes even deeper into intermittent failure patterns.

Secret mission

Build a Flaky Station Detector

Your SPC charts and drift detection caught Station 3 and Station 2. But Station 5 is still slipping through the cracks. Its yield looks fine in aggregate, yet individual DUTs are experiencing random failures that mask a real problem. Can you build a detector that catches intermittently failing stations?

Before you go

Clean Up Your Resources

This project runs entirely locally using Python and Jupyter Notebook, so there are no ongoing cloud costs. Decide whether to keep your resources running, pause them to come back later, or delete them entirely.

Resources you used:

• Jupyter notebook (station_health_monitor.ipynb).
• Installed Python packages: scikit-learn, scipy, pycontrolcharts, plotly.
• Any virtual environment created for this project.

Keep everything running

No action needed. Choose this if you are still actively building or plan to extend the notebook with additional detection modules.

• Your notebook and installed packages will remain available for future use.
• Since everything runs locally, there are no costs associated with keeping your files.

Pause - I'll come back to this later

Shut down running processes to free up memory while keeping your project files intact.

• Stop the Jupyter Notebook server by pressing Ctrl+C in the terminal where it is running.
• Your notebook file and installed packages will still be there when you come back.
• To resume, navigate to your project folder and run jupyter notebook again.

Delete - I don't want to use this again

Remove all project resources and start fresh.

• Stop the Jupyter Notebook server by pressing Ctrl+C in the terminal.
• Delete the notebook file station_health_monitor.ipynb and any associated project folder.
• If you created a virtual environment for this project, delete the virtual environment folder.
• If you installed the packages globally and want to remove them, run pip uninstall scikit-learn scipy pycontrolcharts plotly in your terminal.

Mission accomplished

That's a Wrap!

Nice work! You just built a Test Station Health Monitor that catches what traditional pass/fail testing completely misses. Your notebook can now detect drift, flakiness, and silent yield erosion across a production test floor before any of them become field failures.

Here's what you accomplished:

• Built a synthetic data generator modeling NI TestStand-style output with three hidden failure modes injected across 6 test stations.
• Applied a layered detection pipeline using SPC control charts with Nelson rules, KS drift tests against a golden baseline, and Isolation Forest anomaly scoring to catch multivariate patterns in "passing" DUTs.
• Combined all detection signals into a unified Station Health Report Card with interactive Plotly visualizations that tell a clear, actionable story about test floor health.
• Built a Flaky Station Detector using pass/fail consistency scoring and measurement variance analysis to catch intermittently failing stations that game simple yield metrics.

Ready to quiz yourself?