process-capability-calculator

process-capability-calculator

You are process-capability-calculator - a specialized skill for analyzing process capability with respect to specifications.

Overview

This skill enables AI-powered capability analysis including:

• Capability index calculation (Cp, Cpk)
• Performance index calculation (Pp, Ppk)
• Specification limit analysis
• Normality testing (Shapiro-Wilk, Anderson-Darling)
• Non-normal capability analysis (Box-Cox transformation)
• PPM defect rate estimation
• Capability histogram with distribution overlay
• Six Sigma level calculation

Prerequisites

• Python 3.8+ with numpy, scipy, statsmodels
• Process measurement data
• Specification limits (USL, LSL)

Capabilities

1. Capability Index Calculation (Cp, Cpk)

import numpy as np
from scipy import stats

def calculate_capability_indices(data, usl, lsl, subgroup_size=None):
    """
    Calculate Cp, Cpk capability indices

    Uses within-subgroup variation (R-bar/d2 or S-bar/c4)
    Requires stable process
    """
    x = np.array(data)

    # Process statistics
    x_bar = np.mean(x)
    specification_width = usl - lsl

    # Estimate sigma within (short-term variation)
    if subgroup_size and subgroup_size > 1:
        # Use pooled standard deviation or R-bar method
        # Simplified: using overall std with bias correction
        d2 = {2: 1.128, 3: 1.693, 4: 2.059, 5: 2.326}
        # In practice, calculate R-bar from subgroups
        sigma_within = np.std(x, ddof=1)  # Simplified
    else:
        # For individuals, use moving range
        mr = np.abs(np.diff(x))
        mr_bar = np.mean(mr)
        sigma_within = mr_bar / 1.128  # d2 for n=2

    # Capability indices
    cp = specification_width / (6 * sigma_within)

    cpu = (usl - x_bar) / (3 * sigma_within)
    cpl = (x_bar - lsl) / (3 * sigma_within)
    cpk = min(cpu, cpl)

    return {
        "Cp": round(cp, 3),
        "Cpk": round(cpk, 3),
        "Cpu": round(cpu, 3),
        "Cpl": round(cpl, 3),
        "sigma_within": round(sigma_within, 4),
        "process_mean": round(x_bar, 4),
        "usl": usl,
        "lsl": lsl,
        "interpretation": interpret_capability(cpk)
    }

def interpret_capability(cpk):
    if cpk >= 2.0:
        return "World class (6 sigma)"
    elif cpk >= 1.67:
        return "Excellent (5 sigma)"
    elif cpk >= 1.33:
        return "Good (4 sigma)"
    elif cpk >= 1.0:
        return "Capable (3 sigma)"
    elif cpk >= 0.67:
        return "Poor (2 sigma)"
    else:
        return "Incapable (< 2 sigma)"

2. Performance Index Calculation (Pp, Ppk)

def calculate_performance_indices(data, usl, lsl):
    """
    Calculate Pp, Ppk performance indices

    Uses overall variation (long-term)
    Does not require stable process
    """
    x = np.array(data)

    x_bar = np.mean(x)
    sigma_overall = np.std(x, ddof=1)  # Sample standard deviation
    specification_width = usl - lsl

    # Performance indices
    pp = specification_width / (6 * sigma_overall)

    ppu = (usl - x_bar) / (3 * sigma_overall)
    ppl = (x_bar - lsl) / (3 * sigma_overall)
    ppk = min(ppu, ppl)

    return {
        "Pp": round(pp, 3),
        "Ppk": round(ppk, 3),
        "Ppu": round(ppu, 3),
        "Ppl": round(ppl, 3),
        "sigma_overall": round(sigma_overall, 4),
        "process_mean": round(x_bar, 4),
        "interpretation": interpret_capability(ppk)
    }

3. Normality Testing

def test_normality(data):
    """
    Test data for normality using multiple tests
    """
    x = np.array(data)
    n = len(x)

    results = {}

    # Shapiro-Wilk test (best for n < 50)
    if n <= 5000:
        stat, p_value = stats.shapiro(x)
        results["shapiro_wilk"] = {
            "statistic": round(stat, 4),
            "p_value": round(p_value, 4),
            "conclusion": "Normal" if p_value > 0.05 else "Non-normal"
        }

    # Anderson-Darling test
    ad_result = stats.anderson(x, dist='norm')
    results["anderson_darling"] = {
        "statistic": round(ad_result.statistic, 4),
        "critical_values": dict(zip(
            ['15%', '10%', '5%', '2.5%', '1%'],
            [round(cv, 4) for cv in ad_result.critical_values]
        )),
        "conclusion": "Normal" if ad_result.statistic < ad_result.critical_values[2] else "Non-normal"
    }

    # D'Agostino-Pearson test (n > 20)
    if n >= 20:
        stat, p_value = stats.normaltest(x)
        results["dagostino_pearson"] = {
            "statistic": round(stat, 4),
            "p_value": round(p_value, 4),
            "conclusion": "Normal" if p_value > 0.05 else "Non-normal"
        }

    # Overall assessment
    normal_count = sum(1 for r in results.values() if r.get("conclusion") == "Normal")
    results["overall_assessment"] = "Normal" if normal_count >= 2 else "Non-normal"

    # Descriptive statistics
    results["descriptive"] = {
        "skewness": round(stats.skew(x), 3),
        "kurtosis": round(stats.kurtosis(x), 3),
        "n": n
    }

    return results

4. Non-Normal Capability Analysis

from scipy.stats import boxcox
from scipy.optimize import brentq

def nonnormal_capability(data, usl, lsl, method='percentile'):
    """
    Calculate capability for non-normal data

    Methods:
    - percentile: Use empirical percentiles
    - boxcox: Transform to normal using Box-Cox
    - weibull: Fit Weibull distribution
    """
    x = np.array(data)

    if method == 'percentile':
        # Percentile method (distribution-free)
        p0135 = np.percentile(x, 0.135)  # Lower 0.135%
        p99865 = np.percentile(x, 99.865)  # Upper 99.865%
        median = np.median(x)

        # Equivalent Cp
        spread = p99865 - p0135
        cp_equiv = (usl - lsl) / spread if spread > 0 else np.inf

        # Equivalent Cpk
        cpu_equiv = (usl - median) / ((p99865 - median) * 2) if (p99865 - median) > 0 else np.inf
        cpl_equiv = (median - lsl) / ((median - p0135) * 2) if (median - p0135) > 0 else np.inf
        cpk_equiv = min(cpu_equiv, cpl_equiv)

        return {
            "method": "percentile",
            "Cp_equivalent": round(cp_equiv, 3),
            "Cpk_equivalent": round(cpk_equiv, 3),
            "p0135": round(p0135, 4),
            "p99865": round(p99865, 4),
            "median": round(median, 4)
        }

    elif method == 'boxcox':
        # Box-Cox transformation
        # Shift data if necessary (must be positive)
        shift = 0
        if np.min(x) <= 0:
            shift = abs(np.min(x)) + 1
            x_shifted = x + shift
        else:
            x_shifted = x

        # Find optimal lambda
        transformed, lambda_opt = boxcox(x_shifted)

        # Transform specification limits
        if lambda_opt == 0:
            usl_t = np.log(usl + shift)
            lsl_t = np.log(lsl + shift)
        else:
            usl_t = ((usl + shift)**lambda_opt - 1) / lambda_opt
            lsl_t = ((lsl + shift)**lambda_opt - 1) / lambda_opt

        # Calculate capability on transformed data
        result = calculate_capability_indices(transformed, usl_t, lsl_t)
        result["method"] = "boxcox"
        result["lambda"] = round(lambda_opt, 4)
        result["shift"] = shift

        return result

    return None

5. PPM and Sigma Level Calculation

def calculate_ppm_sigma(cpk, process_centered=True):
    """
    Calculate expected PPM defect rate and sigma level
    """
    if process_centered:
        # Symmetric distribution around target
        z = 3 * cpk
        ppm_total = 2 * (1 - stats.norm.cdf(z)) * 1e6
        ppm_upper = ppm_total / 2
        ppm_lower = ppm_total / 2
    else:
        # Use Cpk for worst side
        z = 3 * cpk
        ppm_worst = (1 - stats.norm.cdf(z)) * 1e6
        ppm_total = ppm_worst  # Approximate

    # Sigma level (with 1.5 sigma shift convention)
    sigma_short_term = 3 * cpk
    sigma_long_term = sigma_short_term + 1.5  # Industry convention

    return