Reasoning Under Uncertainty · Probabilistic methods for condition monitoring

Interactive explorers

Three live tools, one per theme. Tune a kernel bandwidth and watch a density over- or under-smooth, accumulate evidence and watch a fault posterior climb, or sweep the Weibull shape and watch the failure regime change.

Distribution and anomaly, what does normal look like, and what is an outlier?

Before any diagnosis, condition monitoring needs a model of normal: where the healthy data live across operating regimes, so a new point can be scored as familiar or novel, and a single reading can be judged an outlier without distorting the very statistics meant to catch it.

Kernel density estimation (KDE)

distribution

A nonparametric estimate of the data's probability density: place a smooth bump on every observation and sum them. It needs no distributional assumption and yields a novelty score, the negative log density, that flags points unlike anything seen before.

In plain words

A histogram answers where the data are, but crudely, with arbitrary bins and hard edges. Kernel density estimation smooths that into a continuous density by centering a small Gaussian on each observation and adding them up, so dense regions rise and empty regions stay near zero. The one real choice is the bandwidth, the width of those bumps: too narrow and the estimate is spiky and chases noise, too wide and real structure like separate operating regimes is blurred away. Once fitted, the density turns into an anomaly detector for free, since the negative log density is high exactly where a new point is unlike the training data, no fault labels required.

Health-management example

A pump's healthy behavior occupies two clusters in the load-temperature plane, idle and full load. A kernel density estimate over historical healthy data captures both; a new reading at high temperature but low load falls in a near-empty region, earning a high novelty score and flagging an abnormal operating state.

import numpy as np
from scipy import stats

def kde_pdf(X_train, X_eval, bw=None):
    kde = stats.gaussian_kde(np.atleast_2d(X_train.T), bw_method=bw)
    return kde(np.atleast_2d(X_eval.T))

def novelty_score(X_train, X_eval, bw=None):
    return -np.log(kde_pdf(X_train, X_eval, bw) + 1e-300)

Practical note

Everything hinges on the bandwidth; cross-validate it rather than trusting a default, and remember the rule-of-thumb bandwidths assume roughly Gaussian data. KDE degrades in high dimensions, where data become sparse and every point looks novel, so reduce dimensionality or model per-regime. Standardize features first, since one wide-ranging variable will otherwise dominate the distance.

Silverman (1986), Density Estimation; Scott (2015), Multivariate Density Estimation.

Gaussian mixture models (GMM)

distribution

A density modelled as a weighted blend of Gaussian components, fitted by expectation-maximization. It captures several distinct operating regimes at once, assigns each point a soft membership, and gives a parametric likelihood for anomaly scoring.

In plain words

When healthy operation has several modes, a single Gaussian is wrong and a kernel estimate, while flexible, gives no labels. A Gaussian mixture splits the difference: it models the density as a few Gaussian components, each a regime with its own mean, shape, and weight, and fits them by expectation-maximization, which alternates between softly assigning points to components and re-estimating each component from its assigned points. The result both clusters the data into regimes and provides a smooth likelihood, so a new point gets a regime membership and a likelihood, with low likelihood marking an anomaly. The number of components is chosen by an information criterion.

Health-management example

A turbine runs in three regimes, startup, steady, and peaking. A three-component mixture recovers them as separate Gaussians, so each new operating point is assigned to its regime and scored; a point that fits no component well, low under all three, is flagged as an off-nominal state worth investigating.

# see scripts/probstats.py: gmm_fit(X, k) runs EM
model = gmm_fit(X, k=3)            # weights, means, covs, labels
scores = gmm_score(model, X_new)  # per-point log-likelihood
anomaly = scores < threshold      # low likelihood = off-nominal

Practical note

Choose the number of components by BIC rather than by eye, and regularize the covariances or expectation-maximization can collapse a component onto a single point with infinite likelihood. The fit is only locally optimal and depends on initialization, so restart several times and keep the best. Full covariances are expressive but data-hungry; in high dimensions prefer diagonal or shared covariances.

Bishop (2006), Pattern Recognition and ML, ch. 9; McLachlan & Peel (2000), Finite Mixture Models.

Robust outlier tests

distribution

Tests that flag outliers without being corrupted by them. The MAD-based z-score replaces the mean and standard deviation with the median and median absolute deviation; Grubbs and the generalized ESD test give significance for one or several outliers in a near-normal sample.

In plain words

The obvious outlier rule, more than three standard deviations from the mean, has a fatal flaw: a real outlier inflates both the mean and the standard deviation it is measured against, so it hides itself, a problem called masking. Robust statistics fix this by using estimators outliers cannot easily distort. The modified z-score uses the median and the median absolute deviation, which a few extreme points barely move, and flags scores beyond about 3.5. When a formal test is wanted, Grubbs' test judges the single most extreme point against a t-distribution critical value, and the generalized extreme studentized deviate test extends this to several outliers, removing and re-testing so that multiple outliers do not mask one another.

Health-management example

A temperature channel has two stuck-high spikes and one dropout. The classical three-sigma band, widened by the spikes, misses the dropout; the MAD-based score, anchored on the median, flags all three, and the generalized ESD test confirms three significant outliers rather than one.

import numpy as np

def mad_zscore(x):
    med = np.median(x); mad = np.median(np.abs(x - med))
    z = 0.6745 * (x - med) / (mad + 1e-12)
    return z, np.abs(z) > 3.5      # robust to masking
# grubbs(x), generalized_esd(x, max_out) in scripts/probstats.py

Practical note

Grubbs and ESD assume the bulk of the data is roughly normal, so check that before trusting their p-values; for skewed data prefer the distribution-free MAD score or quantile rules. Set the maximum number of outliers in ESD generously, since under-setting it re-introduces masking. Treat a flagged outlier as a hypothesis, not a verdict; it may be a real transient worth keeping, not noise to discard.

Rousseeuw & Hubert (2011), WIREs Data Mining; Rosner (1983), Technometrics 25.

Probabilistic reasoning, what does the evidence imply about a hidden fault?

A fault is hidden; only noisy symptoms are seen. Probabilistic inference runs the logic backward from evidence to cause, carries belief forward as data accumulate, and estimates a hidden state from a stream of measurements, the input to model-based residual generation.

Bayesian updating

probabilistic

Bayes' rule as the engine of diagnosis: combine a prior belief with the likelihood of the evidence to get a posterior, and repeat as evidence arrives. Belief in each fault hypothesis sharpens with every observation.

In plain words

Diagnosis is inference from effect to cause, exactly what Bayes' rule formalizes. It multiplies the prior probability of each hypothesis by how likely the observed evidence would be if that hypothesis were true, then renormalizes, giving the posterior probability of each cause given the evidence. Applied repeatedly, with each posterior becoming the next prior, it accumulates evidence over time, so a weak prior plus a run of consistent observations compounds into near-certainty while a single ambiguous reading shifts belief only a little. For conjugate cases like estimating a defect rate, the update is a closed-form Beta-Binomial step with a credible interval, no simulation needed.

Health-management example

A bearing is healthy with prior probability 0.9. Each vibration snapshot is more consistent with an outer-race fault than with health; after a dozen such snapshots the posterior for the outer-race hypothesis rises from 0.1 toward 0.97, isolating the fault from its competitors as evidence accumulates.

import numpy as np

def bayes_update(prior, likelihood):
    post = np.asarray(prior) * np.asarray(likelihood)
    return post / post.sum()

post = prior
for lik in evidence_stream:        # accumulate evidence
    post = bayes_update(post, lik)

Practical note

The prior matters when evidence is scarce, so set it from reliability data or fleet base rates rather than reflexively choosing uniform; state it explicitly so others can challenge it. The likelihoods must be calibrated, since a mis-specified likelihood drives the posterior confidently wrong. Beware treating correlated observations as independent, which over-counts evidence and produces false certainty.

Gelman et al. (2013), Bayesian Data Analysis; Jaynes (2003), Probability Theory.

Bayesian networks

probabilistic

A graph whose nodes are faults and symptoms and whose edges encode causal and conditional-independence structure. Given observed symptoms it infers the posterior over hidden faults, isolating the likely cause, the probabilistic cousin of partial-correlation reasoning.

In plain words

When many faults and symptoms interact, flat Bayesian updating becomes unwieldy and ignores structure. A Bayesian network draws that structure as a directed graph: each fault points to the symptoms it causes, and the graph's missing edges assert conditional independencies that make the joint distribution compact and inference tractable. Each node carries a conditional probability table, the chance it is present given its parents. Clamping the observed symptoms and inferring the unobserved fault nodes yields the posterior probability of each fault, so the network performs fault isolation: it weighs competing causes against the full pattern of evidence, naturally explaining away one fault when another better accounts for the symptoms.

Health-management example

A network links a bearing fault to both high vibration and high temperature. Observing high vibration alone raises the fault probability from 0.05 to 0.31; adding high temperature, which the same fault explains, pushes it to 0.68, while a competing cause that would explain only temperature is down-weighted by the joint evidence.

# nodes: name -> {parents, cpt: parent-tuple -> P(node=1)}
net = {
  "Fault": {"parents": [], "cpt": {(): 0.05}},
  "Vib":   {"parents": ["Fault"], "cpt": {(0,): 0.1, (1,): 0.85}},
  "Temp":  {"parents": ["Fault"], "cpt": {(0,): 0.15, (1,): 0.7}},
}
p = bayes_net_infer(net, {"Vib": 1, "Temp": 1}, "Fault")  # -> 0.68

Practical note

Inference is exact only for small or sparse networks; dense graphs need approximate methods, and exact enumeration is exponential in the unobserved nodes. The conditional probability tables must come from data or elicited expertise, and a confidently wrong table corrupts the posterior. The structure encodes causal assumptions; validate them, since a missing edge asserts an independence that may not hold.

Pearl (1988), Probabilistic Reasoning; Koller & Friedman (2009), Probabilistic Graphical Models.

Kalman filtering

probabilistic

Recursive Bayesian state estimation for linear-Gaussian systems. It fuses a model-based prediction with each noisy measurement, weighting by uncertainty, to track a hidden state. Its one-step prediction error, the innovation, is the whitened residual that fault detection monitors.

In plain words

Many health variables are hidden states observed only through noise: a true vibration level, a degradation index, a temperature behind a lagging sensor. The Kalman filter estimates them recursively in two steps. It predicts the next state from a dynamics model, inflating its uncertainty, then corrects that prediction with the new measurement, trusting the measurement more when it is precise and the model more when it is not, an optimal blend for linear-Gaussian systems. The byproduct that matters for monitoring is the innovation, the difference between the measurement and its prediction: on a healthy system it is zero-mean white noise, so a persistent or growing innovation is a model-based residual signalling that the system has departed from its assumed behavior.

Health-management example

A degradation index is measured noisily each cycle. The Kalman filter smooths the noisy readings into a clean estimate of the underlying level and its rate of change, and its innovation sequence stays flat while the system tracks the model, then drifts when degradation accelerates beyond the assumed dynamics, raising a residual alarm.

# see scripts/probstats.py: kalman_filter(z, F, H, Q, R, x0, P0)
kf = kalman_filter(z, F, H, Q, R, x0, P0)
state = kf["x_filt"]          # filtered hidden state
resid = kf["innov"]           # innovation = whitened residual

Practical note

The filter is only as good as its models; tune the process and measurement covariances, since too small a process covariance makes it ignore real change and too large makes it chase noise. It is optimal only for linear-Gaussian systems, so use the extended or unscented variants for mild nonlinearity and a particle filter for strong nonlinearity. Check that the innovations are actually white; colored innovations mean the model is mis-specified.

Kalman (1960), J. Basic Eng. 82; Simon (2006), Optimal State Estimation.

Particle filtering

probabilistic

Recursive Bayesian estimation for nonlinear, non-Gaussian systems. It represents the belief by a swarm of weighted samples, propagates them through the true dynamics, weights them by each measurement, and resamples, tracking a hidden health state where the Kalman filter's assumptions fail.

In plain words

Degradation is often nonlinear, a crack that grows faster as it lengthens, a wear process that accelerates, and the Kalman filter's linear-Gaussian assumptions break down. A particle filter drops those assumptions by representing the belief about the state not as a single Gaussian but as a cloud of weighted samples, the particles. Each step it pushes every particle through the actual nonlinear dynamics, reweights them by how well they predict the new measurement, and resamples to concentrate particles where the belief is high. The weighted cloud approximates the full posterior, so it captures skew and multiple modes, and its spread is an honest uncertainty band, valuable for projecting a degradation state and its confidence forward to a remaining-useful-life estimate.

Health-management example

A fatigue crack grows nonlinearly and is measured indirectly and noisily. A particle filter tracks the hidden crack length, its particle cloud widening as uncertainty grows between inspections and tightening at each measurement, and the 5th-to-95th-percentile band projected forward gives a probabilistic remaining-useful-life estimate.

# see scripts/probstats.py: particle_filter(z, f_trans, h_obs, ...)
pf = particle_filter(z, f_trans, h_obs, q_std, r_std, x0)
est = pf["mean"]              # posterior mean health state
band = (pf["lo"], pf["hi"])   # 5-95% band -> RUL uncertainty

Practical note

Particle filters need enough particles, and the count must grow with state dimension or the filter degenerates, the curse of dimensionality. Watch for weight collapse, where one particle takes all the weight; resampling and a small process-noise jitter mitigate it. The dynamics and measurement models still must be right, and the estimate is stochastic, so fix the seed for reproducibility.

Gordon, Salmond & Smith (1993), IEE Proc. F 140; Arulampalam et al. (2002), IEEE TSP 50.

Reliability and prognostics, how long until failure, and at what risk?

Detection and state estimation say what is wrong now; reliability analysis says what it means for the future. From failure times, with censoring handled honestly, it estimates the time-to-failure distribution, the survival curve, the effect of operating conditions, and the current failure rate.

Weibull analysis

reliability

The workhorse life distribution. Its shape parameter names the failure regime, below one infant mortality with a falling hazard, one random failure, above one wear-out with a rising hazard, and its scale is the characteristic life. A probability plot fits and checks it at once.

In plain words

The Weibull distribution is the default model for time-to-failure because its two parameters carry physical meaning. The shape parameter distinguishes the three classic regimes: less than one means the hazard, the failure rate among survivors, falls over time, the signature of infant mortality from manufacturing defects; equal to one means a constant hazard, random failure independent of age; greater than one means a rising hazard, wear-out from fatigue or erosion. The scale parameter is the characteristic life, the age by which about 63 percent have failed. Fitting is by maximum likelihood, and a Weibull probability plot, which linearizes the distribution, doubles as a goodness-of-fit check: the points fall on a line when the model fits.

Health-management example

A fleet of bearings is fitted with a Weibull shape of about 2.3, confirming a wear-out mechanism rather than random failure. The reliability curve then gives the probability a unit survives to any age, and the B10 life, the age at which ten percent have failed, sets a conservative maintenance interval.

import numpy as np
from scipy import stats

def weibull_fit(t):
    shape, loc, scale = stats.weibull_min.fit(t, floc=0)
    return dict(shape=shape, scale=scale)          # k<1,=1,>1 regimes

def weibull_reliability(t, k, lam):
    return dict(R=np.exp(-(t/lam)**k), h=(k/lam)*(t/lam)**(k-1))

Practical note

Estimate parameters by maximum likelihood and report confidence intervals, since small failure samples give wide uncertainty. A single Weibull assumes one failure mode; mixed modes, a falling then rising hazard, need a competing-risks or mixture model, visible as a kink in the probability plot. Always include censored units in the fit, since dropping survivors biases the life estimate short.

Weibull (1951), J. Appl. Mech. 18; Meeker & Escobar (1998), Statistical Methods for Reliability.

Kaplan-Meier estimator

reliability

The nonparametric survival curve. It estimates the probability of surviving past each time without assuming any distribution, and it correctly handles censoring, units still running or withdrawn before failing, which naive averages of failure times get wrong.

In plain words

Sometimes you do not want to assume a Weibull or any other form; you just want the empirical probability of survival over time. The Kaplan-Meier estimator provides it as a product over observed failure times, at each one multiplying the running survival by the fraction of the at-risk population that survived that instant. Its essential feature is the correct treatment of censoring: units still operating, or removed for unrelated reasons, contribute to the at-risk count up to their last-seen time and then drop out without counting as failures, so the curve uses all the information without the bias that ignoring or mishandling survivors would introduce. Comparing two such curves, for two duty cycles or designs, shows directly which survives longer.

Health-management example

Two pump populations run under mild and harsh duty, many still operating at study end. Kaplan-Meier curves, with the censored survivors marked, show the harsh-duty population's survival dropping below 50 percent far earlier, quantifying the duty effect without assuming a failure distribution.

import numpy as np

def kaplan_meier(durations, events):   # events: 1 fail, 0 censored
    S = 1.0; t_out = [0.0]; S_out = [1.0]
    for tt in np.unique(durations[events == 1]):
        at_risk = np.sum(durations >= tt)
        d = np.sum((durations == tt) & (events == 1))
        S *= (1 - d / at_risk); t_out.append(tt); S_out.append(S)
    return np.array(t_out), np.array(S_out)

Practical note

Censoring is assumed non-informative; if units are pulled precisely because they look about to fail, the curve is biased, so check why censoring occurred. The estimate is uncertain where few remain at risk, in the long right tail, so show confidence bands and read the tail cautiously. To test whether two curves truly differ, use a log-rank test rather than eyeballing the gap.

Kaplan & Meier (1958), JASA 53; Klein & Moeschberger (2003), Survival Analysis.

Cox proportional hazards

reliability

A regression for failure risk. It relates covariates, operating conditions or sensor features, to the hazard without specifying a baseline distribution, returning a hazard ratio per covariate: the factor by which it multiplies the risk of failure.

In plain words

Reliability rarely depends on age alone; load, temperature, and condition indicators all shift the risk. The Cox proportional-hazards model captures that without committing to a life distribution. It assumes each covariate multiplies an unspecified baseline hazard by a constant factor, so the baseline shape cancels out of the estimation and the covariate effects are fitted from the ordering of failures alone, by maximizing a partial likelihood. Each coefficient becomes a hazard ratio: greater than one means the covariate raises failure risk, less than one means it protects, and one means no effect. It is the standard tool for asking which operating conditions or features actually drive failure, and by how much.

Health-management example

Fitting failure times against vibration, temperature, and load gives hazard ratios of about 2.8 for vibration, 1.7 for temperature, and near one for load. The model thus identifies vibration as the dominant risk driver, each unit increase nearly tripling the hazard, while load shows no significant effect, focusing maintenance attention.

# see scripts/probstats.py: cox_ph(X, durations, events)
fit = cox_ph(X, durations, events)
fit["hazard_ratio"]    # exp(coef): >1 raises risk, <1 protects
fit["coef"], fit["se"] # for confidence intervals on the ratios

Practical note

The core assumption is proportional hazards, that a covariate's effect is constant over time; check it, since a crossing effect violates it and needs a time-varying term or stratification. Hazard ratios describe relative risk, not absolute timing, and the model assumes a log-linear covariate effect, so consider transforms for strongly nonlinear drivers. With many covariates and few failures, regularize to avoid overfitting.

Cox (1972), J. R. Stat. Soc. B 34; Therneau & Grambsch (2000), Modeling Survival Data.

Hazard-rate estimation

reliability

The instantaneous failure rate of a unit still working, and its running total, the cumulative hazard. The Nelson-Aalen estimator gives it nonparametrically; its shape, falling, flat, or rising, locates a population on the bathtub curve.

In plain words

The hazard rate is the chance a surviving unit fails in the next instant, the most diagnostic reliability quantity because its shape reveals the failure mechanism. Plotted over a product's life it traces the bathtub curve: a high but falling hazard early on from infant mortality, a low flat middle of random failures, and a rising tail as wear-out sets in. The Nelson-Aalen estimator builds the cumulative hazard nonparametrically by summing the instantaneous failure fractions at each event time, handling censoring like Kaplan-Meier, and the slope of that cumulative curve is the hazard rate. A steepening curve is the clearest sign that a population has entered wear-out and that preventive replacement is due.

Health-management example

The cumulative hazard of an aging fleet is nearly flat early, then bends sharply upward after a characteristic age, the population entering wear-out. That steepening, rather than any single failure, is what justifies moving from run-to-failure to scheduled replacement before the rising hazard makes failures common.

import numpy as np

def nelson_aalen(durations, events):    # cumulative hazard H(t)
    H = 0.0; t_out = [0.0]; H_out = [0.0]
    for tt in np.unique(durations[events == 1]):
        at_risk = np.sum(durations >= tt)
        d = np.sum((durations == tt) & (events == 1))
        H += d / at_risk; t_out.append(tt); H_out.append(H)
    return np.array(t_out), np.array(H_out)   # slope = hazard rate

Practical note

Nonparametric hazard estimates are noisy, especially the instantaneous rate, so smooth or work with the cumulative hazard whose slope is steadier. The bathtub is a composite of failure modes, not one process, so a real population may show only part of it. Distinguish a rising hazard, genuine wear-out, from a rising failure count driven merely by more units reaching old age; the hazard normalizes by the at-risk population.

Nelson (1972), Technometrics 14; Aalen (1978), Ann. Statist. 6.

Decision map

Start from the question you actually have. The badge marks the family: distribution, probabilistic, or reliability.

If you are asking...	Reach for	Family
Where does normal operation live, across regimes?	KDE / GMM density	distribution
Is this point novel or anomalous?	novelty score (neg log density)	distribution
How many operating regimes are there?	Gaussian mixture	distribution
Is this single reading an outlier?	Grubbs / MAD z-score	distribution
Are several outliers masking each other?	generalized ESD	distribution
What does this evidence imply about a fault?	Bayesian updating	probabilistic
Which fault, given several symptoms?	Bayesian network	probabilistic
What is the hidden state from noisy linear data?	Kalman filter	probabilistic
What is the hidden state under nonlinear dynamics?	particle filter	probabilistic
What is the time-to-failure distribution?	Weibull analysis	reliability
Survival with censored data, no distribution?	Kaplan-Meier	reliability
How do operating conditions affect failure risk?	Cox proportional hazards	reliability
Is the failure rate rising (wear-out)?	hazard / Nelson-Aalen	reliability

Choosing a probabilistic tool

A working order for probabilistic prognostics

Model normal: fit a density or mixture over healthy, regime-aware data and standardize the features.
Score novelty and screen single readings with a robust outlier test; calibrate the threshold conformally.
When a fault is suspected, fuse the symptoms with Bayesian updating or a Bayesian network to isolate the cause.
Estimate the hidden degradation state with a Kalman, unscented, or particle filter, and watch its residual.
Fit a life model (Weibull, or Kaplan-Meier if you avoid a distribution) and read the current hazard.
Project the state to a remaining life with an interval, check calibration, and choose a cost-optimal alarm.

# common stack
pip install numpy scipy
# kde_pdf, gmm_fit, mad_zscore, bayes_net_infer, kalman_filter, particle_filter, weibull_fit, cox_ph, rul_predict, hmm_fit, gp_regression, bocpd, calibration_curve ship in the kit engine

Rules of thumb

Heuristics that hold across most probabilistic condition monitoring, before the method-specific detail.

Model normal before hunting anomalies. A density over healthy, regime-aware data turns novelty into a score; without it, every unusual but benign operating point looks like a fault.
Standardize features first. Density estimates, mixtures, and distance-based scores are all distorted when one variable's range dwarfs another's.
Use robust statistics for outliers. The mean and standard deviation are inflated by the very outliers they should catch; the median and MAD are not.
State the prior explicitly. In Bayesian reasoning the prior matters most when evidence is scarce; draw it from base rates and let others challenge it.
Do not double-count correlated evidence. Treating dependent observations as independent manufactures false certainty in a posterior.
Match the filter to the dynamics. Kalman for linear-Gaussian, extended or unscented for mild nonlinearity, particle for strong nonlinearity or non-Gaussian noise.
Always handle censoring. Units still running carry information; dropping them biases every life and survival estimate toward early failure.
Read the hazard shape, not just the rate. Falling, flat, or rising places a population on the bathtub curve and tells infant mortality from wear-out.

Pick by scenario

Scenario	Pick & why	Family
Flag an abnormal operating point	KDE novelty scoredensity without fault labels	distribution
Separate multi-modal operating regimes	Gaussian mixturesoft regimes + likelihood	distribution
Catch outliers that hide each other	generalized ESDremoves masking	distribution
Isolate a fault from several symptoms	Bayesian networkweighs joint evidence	probabilistic
Accumulate evidence over time	sequential Bayesian updatingposterior as next prior	probabilistic
Track a hidden level from noisy sensors	Kalman filteroptimal linear-Gaussian estimate	probabilistic
Project a nonlinear degradation state	particle filterfull posterior + band	probabilistic
Model time-to-failure	Weibull analysisshape names the regime	reliability
Survival with many units still running	Kaplan-Meierhandles censoring	reliability
Rank what drives failure risk	Cox proportional hazardshazard ratio per covariate	reliability

At a glance

The methods split by the question: what is normal, what does the evidence imply, and how long until failure.

Method	Family	Answers	Output	Key assumption
KDE	distribution	where data live, novelty	density + novelty score	bandwidth choice
GMM	distribution	operating regimes	soft clusters + likelihood	number of components
Robust tests	distribution	is this an outlier	flag + significance	near-normal bulk
Bayesian updating	probabilistic	fault given evidence	posterior probabilities	calibrated likelihoods
Bayesian network	probabilistic	which fault isolates	posterior over faults	graph + CPTs
Kalman filter	probabilistic	hidden state, residual	state + innovation	linear-Gaussian
Particle filter	probabilistic	nonlinear hidden state	state + uncertainty band	enough particles
Weibull	reliability	time-to-failure regime	shape, scale, R(t), h(t)	single failure mode
Kaplan-Meier	reliability	survival over time	survival curve	non-informative censoring
Cox PH	reliability	what drives risk	hazard ratio per covariate	proportional hazards
Hazard / N-A	reliability	is the rate rising	cumulative hazard	non-informative censoring

Worked case study: from anomaly to remaining life on a pump bearing

One unit in a fleet of 240 drifts off-nominal. The chain below detects it, estimates its hidden degradation state, and turns that into a failure probability and a remaining-life figure. Every number is computed by the engine.

Step 1

Detect the anomaly with a density model

A kernel density estimate over the healthy fleet's load-temperature operating points defines normal. The suspect unit, running hot at moderate load, lands in a near-empty region: its novelty score is 25.5, against a fleet median of 5.9 and a 95th percentile of 7.8. The operating point is unlike anything healthy, without any fault label being needed.

Step 2

Estimate the hidden state with a Kalman filter

The unit's noisy vibration-RMS trend is fed to a two-state Kalman filter tracking level and rate. It estimates the current level at 3.32 and a rising rate of 0.061 per cycle. Projecting that rate to the alarm threshold of 4.5 gives a remaining useful life of about 19 cycles, with the innovation sequence confirming the trend is real rather than noise.

Step 3

Quantify the risk with a life model

A Weibull fit on this bearing type's failure history gives a shape of 2.27, a wear-out regime, and a characteristic life of 528. At the unit's current age of 300, past the B10 life of 196, reliability is 0.76 and the hazard is climbing. Fusing the anomaly and the rising trend by Bayes' rule lifts the fault probability from a prior of 0.05 to 0.4.

Verdict

What the chain delivered

A density model flagged an operating point no threshold would have caught; a Kalman filter turned a noisy trend into a state estimate and a remaining-life figure; and a Weibull model placed the unit in wear-out past its B10 life. Bayes' rule fused the evidence into a single fault probability. The reflex order for prognostics is: model normal and score novelty, estimate the hidden state and project it, then quantify the time-to-failure risk, and fuse the three into one decision. Detection says something changed; the state estimate says how fast; the life model says how much time is left.

Going further

Beyond the core three families, a set that makes remaining-useful-life prediction first-class, sharpens anomaly detection, completes the inference toolkit, and turns probabilities into decisions. Two more explorers first.

Remaining useful life: the prognostic capstone

Detection and state estimation say what is wrong; prognostics says how much time is left. These two make remaining-useful-life prediction first-class, and give the metrics to judge whether a prediction is any good.

▸

Remaining useful life prediction

going further

Project a fitted degradation trend forward to the failure threshold and report the crossing time as a remaining life with a prediction interval. The interval widens the further ahead it reaches, an honest statement of forecast uncertainty.

In plain words

A health indicator that drifts toward a failure threshold contains a forecast: extrapolate its trend and the crossing time is the remaining useful life. The value of doing this probabilistically is the interval. Fitting the trend leaves uncertainty in its slope and intercept, and propagating that uncertainty to the crossing gives a band rather than a point, narrow when many consistent observations pin the trend, wide early on or when the signal is noisy. Feeding a filtered state, from the Kalman or particle card, into this projection rather than the raw signal sharpens both the estimate and the band, and the band is what turns a remaining-life number into a maintenance decision with a known confidence.

Health-management example

A bearing's vibration index has climbed steadily for thirty cycles. Extrapolating the fitted trend to the alarm level gives a remaining life of about five cycles with an interval of four to six, tight enough to schedule the replacement during the next planned downtime rather than reacting to a failure.

# see scripts/probstats.py: rul_predict(t, y, threshold)
rl = rul_predict(t, health_index, threshold=4.5)
rl["rul"]                      # remaining cycles to threshold
rl["rul_lo"], rl["rul_hi"]     # 5-95% prediction interval

Practical note

The interval is only as honest as the trend model; a linear fit to a nonlinear degradation will be confidently wrong near the end, so check the trend shape and consider a nonlinear or filtered projection. Report the interval, never the point alone, and recompute as new data arrive. Remember the threshold itself is uncertain, so a conservative threshold buys margin.

Si et al. (2011), Eur. J. Oper. Res. 213; Lei et al. (2018), Mech. Syst. Signal Process. 104.

▸

Prognostic performance metrics

going further

The metrics that judge a remaining-life predictor: alpha-lambda accuracy, whether each prediction sits in a shrinking cone around the truth; prognostic horizon, how early it enters and stays there; and relative accuracy. Raw error alone does not capture what a prognostic needs.

In plain words

A remaining-life predictor is not judged like a regressor, because what matters is being right with enough lead time and getting more accurate as failure approaches. The standard metrics encode that. Alpha-lambda accuracy draws a cone of plus-or-minus a fraction around the true remaining life and asks whether the prediction stays inside it as the unit ages. The prognostic horizon is the earliest point after which the predictions remain in that cone, the lead time the method actually delivers. Relative accuracy summarizes the average closeness. Together they reward a predictor that commits early and converges, and penalize one that is only right at the very end when it is too late to act.

Health-management example

Two predictors have similar average error, but one enters the twenty-percent accuracy cone fifteen cycles before failure and the other only five. The prognostic horizon makes the first clearly better, since it gives three times the lead time to plan maintenance, a distinction raw error hides.

# see scripts/probstats.py: prognostic_metrics(t, rul_pred, rul_true)
pm = prognostic_metrics(t, rul_pred, rul_true, alpha=0.2)
pm["prognostic_horizon"]   # earliest time predictions stay in the cone
pm["rel_accuracy"]         # average closeness to true RUL

Practical note

Pick the alpha and the metrics with the maintenance lead time in mind; a tight cone is meaningless if the required horizon is short. Evaluate on run-to-failure data with known end of life, and average over many units, since a single trajectory is anecdote. Report the horizon alongside accuracy, since an accurate but late prediction has little operational value.

Saxena et al. (2008), Int. Conf. Prognostics & Health Management; Saxena et al. (2010), PHM Society.

Sharper anomaly detection

The opening cards score novelty and flag outliers; these three sharpen the edges: correlation-aware multivariate distance, a principled tail for thresholds, and a calibrated cutoff with a coverage guarantee.

▸

Mahalanobis distance with robust covariance

going further

A multivariate outlier distance that accounts for the correlation between features, with a robust center and covariance so the outliers cannot inflate the very statistics meant to catch them. It is the multivariate, correlation-aware extension of the robust z-score, and the link to the Hotelling chart.

In plain words

A point can be unremarkable in every single variable yet impossible in combination, sitting far off the correlation the features normally obey. Euclidean distance misses this; the Mahalanobis distance does not, because it measures distance in units of the data's covariance, stretching and rotating space so the normal cloud becomes a sphere and off-correlation points stand out. The catch is the same masking problem as in one dimension: outliers inflate the covariance and hide. A robust estimate of the center and covariance, computed by down-weighting or trimming the most extreme points, fixes that, and the resulting robust distances expose multivariate outliers a classical Mahalanobis distance would absorb. It is the natural companion to the Hotelling chart in your multivariate-monitoring toolkit.

Health-management example

Temperature and load are each within their normal range, but the unit is hot for how lightly it is loaded, a combination off the usual correlation. The robust Mahalanobis distance flags it while a classical distance, its covariance stretched toward the very outlier, rates it nearly normal.

# see scripts/probstats.py: mahalanobis(X, robust=True)
m = mahalanobis(X, robust=True)
m["dist"]        # robust Mahalanobis distance per point
m["outlier"]     # flagged against the chi-square threshold

Practical note

Use a robust covariance (a minimum-covariance-determinant style estimate) whenever the training data may contain outliers, which is almost always. Compare distances to a chi-square quantile with the right degrees of freedom for a threshold. In high dimensions the covariance is hard to estimate; regularize it or reduce dimensionality first.

Rousseeuw & Van Zomeren (1990), JASA 85; Hubert et al. (2018), WIREs Comput. Stat.

▸

Extreme value theory for thresholds

going further

A principled alarm level from the tail of a health index. Exceedances over a high threshold follow a generalized Pareto distribution; fitting it sets the level a healthy signal crosses only with a chosen small probability, replacing an arbitrary number of sigmas.

In plain words

Alarm thresholds are too often set at three sigmas or a round number, neither of which corresponds to a stated false-alarm rate, especially when the health index is not Gaussian. Extreme value theory addresses exactly the tail that matters. Its peaks-over-threshold result says the exceedances above a high threshold converge to a generalized Pareto distribution regardless of the underlying shape, so fitting that distribution to the observed exceedances lets you read off the level that a healthy signal will exceed only with, say, a one-in-a-thousand probability. The alarm is then expressed as a return level with a meaning, rather than a heuristic, and it adapts to heavy tails that a Gaussian assumption would badly underestimate.

Health-management example

A health index has occasional benign spikes, so a three-sigma alarm fires too often. Fitting a generalized Pareto to its tail and choosing a one-in-a-thousand exceedance level sets a threshold that tolerates the benign spikes while still catching genuinely extreme excursions.

# see scripts/probstats.py: pot_threshold(x, q=0.95, exceed_prob=1e-3)
pot = pot_threshold(health_index, q=0.95, exceed_prob=1e-3)
pot["level"]     # alarm level for the chosen exceedance probability
pot["shape"]     # GPD shape: >0 heavy tail, <0 bounded

Practical note

The choice of the threshold for defining exceedances is a bias-variance tradeoff: too low and the asymptotics do not hold, too high and too few exceedances remain to fit; inspect a mean-residual-life plot. The fit needs enough tail data, so it is for well-sampled indices, not rare events with a handful of points. Extrapolating far beyond the data carries real uncertainty, so quote confidence on the return level.

Coles (2001), Introduction to Statistical Modeling of Extreme Values; Pickands (1975), Ann. Statist. 3.

▸

Conformal anomaly detection

going further

A wrapper that turns any anomaly score into a calibrated alarm with a finite-sample coverage guarantee. Using scores on clean calibration data, it sets the cutoff so a normal point exceeds it with probability at most alpha, no distributional assumption required.

In plain words

Any of the scores in this guide, a kernel density novelty, a mixture likelihood, a Mahalanobis distance, leaves one question open: where to put the threshold so the false-alarm rate is actually what you claim. Conformal prediction answers it with a guarantee. Given a set of anomaly scores computed on data known to be normal, it sets the cutoff at the appropriate empirical quantile, with a small finite-sample correction, so that a future normal point exceeds it with probability at most the chosen alpha, exactly, regardless of the score's distribution or the data's. It converts an uncalibrated score into a test with a controlled false-alarm rate, which is what a deployable alarm needs and what most anomaly detectors lack.

Health-management example

A kernel-density novelty score has no natural threshold. Calibrating it conformally on a month of known-healthy data sets a cutoff guaranteed to false-alarm on at most five percent of normal points, so the alarm rate in deployment matches the stated tolerance rather than drifting with the score's scale.

# see scripts/probstats.py: conformal_threshold(cal_scores, alpha=0.05)
ct = conformal_threshold(scores_on_normal_data, alpha=0.05)
alarm = new_score > ct["threshold"]   # <=5% false-alarm guarantee

Practical note

The guarantee assumes the calibration and future data are exchangeable, so it breaks under distribution shift; re-calibrate when the operating regime changes. It controls the marginal false-alarm rate, not the per-point one, and says nothing about detection power, which still depends on the underlying score. Hold out a clean calibration set separate from training.

Vovk et al. (2005), Algorithmic Learning in a Random World; Angelopoulos & Bates (2023), Found. Trends ML.

More state and inference tools

Beyond the linear Kalman filter and the bootstrap particle filter, four members that round out the inference family: discrete degradation stages, derivative-free nonlinear filtering, smooth trends with uncertainty, and online fault-onset timing.

▸

Hidden Markov models

going further

A model of a system passing through hidden discrete stages, healthy, degraded, failing, each emitting noisy observations. Baum-Welch learns the stages and their transition probabilities; Viterbi decodes the most likely stage sequence, segmenting a noisy signal into a clean degradation history.

In plain words

Degradation is often naturally discrete: a machine is healthy, then degraded, then failing, with noisy sensor values that overlap between stages. A hidden Markov model captures this as a Markov chain over hidden stages that you cannot see directly, each emitting observations with its own distribution. The Baum-Welch algorithm, an expectation-maximization, learns the stage means, their spreads, and the probabilities of moving between them, all from unlabelled data. The Viterbi algorithm then finds the single most likely sequence of hidden stages behind the observations, turning a jittery signal into a clean segmentation that says when the system crossed from one stage to the next, the discrete-state counterpart to the continuous filters.

Health-management example

A gearbox health signal wanders noisily upward over months. A three-state hidden Markov model, fitted without labels, recovers a healthy, a degraded, and a failing stage, and Viterbi decoding marks the two transition times, giving a maintenance planner a clean stage history rather than a noisy trace.

# see scripts/probstats.py: hmm_fit(obs, n_states); hmm_viterbi(obs, model)
model = hmm_fit(health_signal, n_states=3)   # states sorted by mean
stages = hmm_viterbi(health_signal, model)   # 0=healthy ... 2=failing

Practical note

Choose the number of stages with care and a physical rationale; too many fragments the history, too few merges distinct regimes. Baum-Welch finds only a local optimum, so initialize sensibly and restart. Label-switching means the learned states are unordered, so sort them by mean before interpreting them as healthy-to-failing.

Rabiner (1989), Proc. IEEE 77; Baum et al. (1970), Ann. Math. Statist. 41.

▸

Unscented Kalman filter

going further

Nonlinear state estimation without Jacobians. It propagates a small set of deterministically chosen sigma points through the true nonlinear dynamics and recovers the mean and covariance, more accurate than linearizing and bridging the linear Kalman filter and the particle filter.

In plain words

Between the linear Kalman filter and the fully general particle filter sits a large class of mildly nonlinear, roughly Gaussian problems. The extended Kalman filter handles them by linearizing the dynamics, which needs Jacobians and degrades when the nonlinearity is sharp. The unscented Kalman filter takes a cleaner route: it picks a handful of sigma points that capture the state's mean and covariance, pushes each through the exact nonlinear function, and reconstructs the transformed mean and covariance from the results. This unscented transform is accurate to higher order than linearization, needs no derivatives, and costs only a few function evaluations, making it the practical default for nonlinear-Gaussian state estimation in diagnostics.

Health-management example

A sensor reads the square root of a degradation state, a nonlinear observation. The unscented Kalman filter tracks the underlying state cleanly through that nonlinearity, where a linearized filter would bias the estimate, and without the derivative the extended filter would require.

# see scripts/probstats.py: ukf_filter(z, f, h, Q, R, x0, P0)
uk = ukf_filter(z, f_dynamics, h_observe, Q, R, x0, P0)
uk["x_filt"]     # nonlinear state estimate
uk["innov"]      # innovation residual

Practical note

It assumes the posterior stays roughly Gaussian, so for strongly multimodal or heavy-tailed problems use a particle filter. The sigma-point spread parameters rarely need tuning but can matter for high dimension. As with any Kalman variant, the process and measurement covariances must be set sensibly and the innovations checked for whiteness.

Julier & Uhlmann (2004), Proc. IEEE 92; Wan & Van Der Merwe (2000), IEEE AS-SPCC.

▸

Gaussian-process regression

going further

A nonparametric trend model that returns a calibrated uncertainty band. It fits a smooth function to a degradation signal and, crucially, widens its confidence where data are sparse and as it extrapolates, exactly the behavior wanted for projecting a health trend toward failure.

In plain words

Fitting a degradation trend with a fixed parametric curve forces a shape and gives no honest uncertainty. A Gaussian process instead places a distribution over smooth functions and conditions it on the data, yielding not just a mean trend but a full predictive distribution at every point. Its defining feature for prognostics is the uncertainty: the band is tight where observations are dense and widens where they are sparse and, most importantly, as the prediction extrapolates beyond the last measurement. That growing band is a faithful statement of how little the data constrain the future, so when the trend is projected toward a failure threshold the crossing time inherits a principled interval rather than a false-confident point.

Health-management example

A degradation signal is fitted with a Gaussian process whose mean tracks the trend and whose two-sigma band hugs the data, then flares once the prediction extends past the last observation. Projecting the band to the failure threshold gives a remaining-life interval that honestly reflects the thinning evidence.

# see scripts/probstats.py: gp_regression(t_train, y_train, t_test, ...)
gp = gp_regression(t, y, t_future, length=2.5, sigma_f=2.0, sigma_n=0.2)
gp["mean"], gp["std"]    # trend and uncertainty (std grows ahead)

Practical note

The kernel and its length scale encode the smoothness assumption; choose them by marginal likelihood, not by eye, and a wrong length scale gives over- or under-confident bands. The standard implementation costs cubically in the number of points, so thin or approximate for long histories. Extrapolation reverts to the prior mean far out, so trust the band, not the mean, in the deep future.

Rasmussen & Williams (2006), Gaussian Processes for Machine Learning; Richardson et al. (2017), Mech. Syst. Signal Process.

▸

Bayesian online changepoint detection

going further

An online estimator of the time since the last regime change. It maintains a posterior over the run length using a conjugate predictive, and the run length collapses to zero exactly when the signal's statistics shift, flagging a fault onset as it happens, the probabilistic timing tool alongside your change-detection methods.

In plain words

Knowing that a fault exists is one thing; knowing precisely when it started is another, and it matters for root cause and for resetting baselines. Bayesian online changepoint detection answers it sequentially. At each step it maintains a probability distribution over the run length, the number of steps since the last changepoint, updating it with a conjugate predictive model of the data. While the regime is stable the most likely run length grows steadily; the instant the underlying mean or variance shifts, the new observation fits the grown hypothesis poorly, the run-length posterior collapses toward zero, and a changepoint is declared. Because it is online and probabilistic, it reports not just that a change occurred but when, with a posterior over the timing.

Health-management example

A vibration level shifts to a new baseline midway through a record. The run-length estimate climbs steadily, then drops sharply to near zero at the shift, pinpointing the onset time so the event can be correlated with a load change or a maintenance action.

# see scripts/probstats.py: bocpd(x, hazard=0.01)
bc = bocpd(signal, hazard=0.01)
# expected run length resets to ~0 at a changepoint
onset = int(bc["exp_run_length"].argmin())

Practical note

The hazard rate encodes how often changepoints are expected; set it from domain knowledge, since too high makes the detector jumpy and too low makes it sluggish. The conjugate model assumes a form for the observations (here Gaussian), so transform or generalize it for other data. There is an inherent detection lag, a few samples, while evidence for the change accumulates.

Adams & MacKay (2007), arXiv:0710.3742; Fearnhead & Liu (2007), J. R. Stat. Soc. B 69.

From probabilities to decisions

A probability is only useful if it is honest and if it drives an action. These three close the loop: checking that a stated probability is calibrated, turning it into a cost-optimal alarm, and fusing several detectors into one belief.

▸

Calibration and proper scoring

going further

The check that a model's stated probabilities mean what they say. A reliability diagram compares predicted probability to observed frequency, and the Brier score summarizes calibration and sharpness. A confident-but-wrong model is worse than an honest one.

In plain words

Every probabilistic method in this guide outputs a number between zero and one, but a number is only trustworthy if it is calibrated: of all the times the model says seventy percent, a fault should actually occur about seventy percent of the time. A reliability diagram tests this directly by binning predictions and plotting the stated probability against the observed frequency, with calibration meaning the points lie on the diagonal. The Brier score, the mean squared difference between probability and outcome, condenses both calibration and sharpness into one proper score, one that is optimized only by reporting the true probability. Checking calibration matters because an overconfident model triggers the wrong maintenance decisions even when its rankings are fine, and recalibration can fix it without retraining.

Health-management example

A fault classifier looks accurate but its reliability diagram bows away from the diagonal: when it says ninety percent, faults occur only seventy percent of the time. Recognizing this overconfidence, and recalibrating, prevents the cost model downstream from over-reacting to inflated probabilities.

# see scripts/probstats.py: calibration_curve(probs, labels)
cc = calibration_curve(fault_probs, outcomes, n_bins=10)
cc["conf"], cc["acc"]    # diagonal = calibrated
cc["brier"]              # proper score: lower is better

Practical note

Calibration and accuracy are different; a model can rank perfectly yet be badly calibrated, and the fix is recalibration, not retraining. Use enough bins to see structure but enough data per bin to be stable, and prefer a proper score like Brier or log-loss over accuracy for probabilistic models. Recalibrate after any distribution shift.

Brier (1950), Mon. Weather Rev. 78; Guo et al. (2017), ICML (modern calibration).

▸

Cost-optimal decision thresholds

going further

Turn a fault probability into a maintenance action by choosing the alarm threshold that minimizes expected cost, given the asymmetric costs of a false alarm and a missed fault. It replaces an arbitrary cutoff with one tuned to the real consequences.

In plain words

A calibrated fault probability still has to become a yes-or-no decision, and the right threshold is rarely one-half. The costs are asymmetric: an unnecessary inspection wastes a little, an unexpected failure can be catastrophic, so the threshold should reflect that ratio. Casting it as expected-cost minimization makes this explicit: for each candidate threshold, count the false alarms and the missed faults the decisions would produce, weight them by their costs, and choose the threshold with the lowest expected cost. The result moves the alarm lower when misses are expensive, catching more faults at the price of more false alarms, and higher when false alarms dominate. It is the bridge from a probability to an operational policy.

Health-management example

With a missed bearing failure costing twenty times an unnecessary inspection, the cost-optimal threshold drops to about a quarter rather than a half, so the system alarms earlier and accepts more false positives, because the rare missed failure is what truly hurts.

# see scripts/probstats.py: cost_optimal_threshold(probs, labels, c_fp, c_fn)
co = cost_optimal_threshold(fault_probs, outcomes, c_fp=1.0, c_fn=20.0)
co["best"]        # cost-minimizing alarm threshold

Practical note

The decision is only as good as the cost numbers; elicit them honestly and do a sensitivity analysis, since the optimal threshold can move sharply with the cost ratio. It assumes the probabilities are calibrated, so check that first. Costs may vary by unit and over time, so revisit the threshold rather than fixing it once.

Elkan (2001), IJCAI (cost-sensitive learning); Hand (2009), Mach. Learn. 77.

▸

Evidence fusion (Dempster-Shafer)

going further

A way to combine several detectors into one belief that represents ignorance explicitly. Each source assigns mass to sets of hypotheses, including the whole set as uncertainty, and Dempster's rule fuses them, concentrating belief where they agree, a richer alternative to averaging probabilities.

In plain words

Real diagnosis draws on several detectors, a vibration anomaly, a temperature trend, a model residual, and they must be combined. Simply averaging their probabilities discards how confident each one is and cannot represent plain ignorance. Dempster-Shafer theory lets each source distribute belief mass not only over individual hypotheses but over sets of them, assigning leftover mass to the whole frame as explicit uncertainty. Dempster's rule of combination then fuses two sources by intersecting their focal sets and renormalizing away the conflict, so agreement sharpens belief while disagreement is surfaced rather than hidden. The result is a fused belief that grows more committed than any single source when they concur, and an honest measure of conflict when they do not, useful when sensors are partially reliable or sometimes silent.

Health-management example

A vibration sensor and a temperature sensor each lean toward a fault but reserve some mass for uncertainty. Combining them by Dempster's rule raises the fused belief in the fault above either alone, while a high conflict score would instead warn that the two sensors disagree and should not be naively merged.

# see scripts/probstats.py: ds_combine(m1, m2)
F, H, U = frozenset(["F"]), frozenset(["H"]), frozenset(["F", "H"])
m1 = {{F: 0.55, H: 0.10, U: 0.35}}   # vibration sensor
m2 = {{F: 0.50, H: 0.15, U: 0.35}}   # temperature sensor
fused = ds_combine(m1, m2)         # fused belief in the fault

Practical note

High conflict between sources makes Dempster's rule behave counterintuitively, the classic Zadeh paradox, so monitor the conflict and consider alternative combination rules when it is large. Eliciting the mass functions is the hard part and encodes each sensor's reliability. When all sources give well-calibrated probabilities, Bayesian model averaging is often the simpler and more defensible choice.

Shafer (1976), A Mathematical Theory of Evidence; Dempster (1967), Ann. Math. Statist. 38.

References

Silverman (1986), Density Estimation; Scott (2015), Multivariate Density Estimation.
Bishop (2006), Pattern Recognition and ML, ch. 9; McLachlan & Peel (2000), Finite Mixture Models.
Rousseeuw & Hubert (2011), WIREs Data Mining; Rosner (1983), Technometrics 25.
Gelman et al. (2013), Bayesian Data Analysis; Jaynes (2003), Probability Theory.
Pearl (1988), Probabilistic Reasoning; Koller & Friedman (2009), Probabilistic Graphical Models.
Kalman (1960), J. Basic Eng. 82; Simon (2006), Optimal State Estimation.
Gordon, Salmond & Smith (1993), IEE Proc. F 140; Arulampalam et al. (2002), IEEE TSP 50.
Weibull (1951), J. Appl. Mech. 18; Meeker & Escobar (1998), Statistical Methods for Reliability.
Kaplan & Meier (1958), JASA 53; Klein & Moeschberger (2003), Survival Analysis.
Cox (1972), J. R. Stat. Soc. B 34; Therneau & Grambsch (2000), Modeling Survival Data.
Nelson (1972), Technometrics 14; Aalen (1978), Ann. Statist. 6.

Take it into your editor

The same content is packaged as a VS Code Copilot skill in the progressive-disclosure layout: a routing SKILL.md, per-topic reference docs, a tested probstats.py engine that implements every method with numpy and scipy only, and Copilot integration files (a reliability-analyst chat mode, a /diagnose prompt, and path-scoped instructions). Drop the folder into your kit, or point Copilot at the skill directly.

↓ probabilistic-methods-kit.zip
SKILL.md · references · probstats.py · Copilot agent + prompt

Link resolves when the zip sits beside this file. Run python scripts/probstats.py for the self-test across all three families.

What is normal?

What does the evidence imply?

How long until failure?

Kernel density and the bandwidth

Bayesian updating

The Weibull shape and the failure regime

Kernel density estimation (KDE)

In plain words

Gaussian mixture models (GMM)

In plain words

Robust outlier tests

In plain words

Bayesian updating

In plain words

Bayesian networks

In plain words

Kalman filtering

In plain words

Particle filtering

In plain words

Weibull analysis

In plain words

Kaplan-Meier estimator

In plain words

Cox proportional hazards

In plain words

Hazard-rate estimation

In plain words

Decision map

Choosing a probabilistic tool

A working order for probabilistic prognostics

Rules of thumb

At a glance

Detect the anomaly with a density model

Estimate the hidden state with a Kalman filter

Quantify the risk with a life model

What the chain delivered

Remaining useful life

Calibration

Remaining useful life prediction

In plain words

Prognostic performance metrics

In plain words

Mahalanobis distance with robust covariance

In plain words

Extreme value theory for thresholds

In plain words

Conformal anomaly detection

In plain words

Hidden Markov models

In plain words

Unscented Kalman filter

In plain words

Gaussian-process regression

In plain words

Bayesian online changepoint detection

In plain words

Calibration and proper scoring

In plain words

Cost-optimal decision thresholds

In plain words

Evidence fusion (Dempster-Shafer)

In plain words

References

Take it into your editor