Degradation Modelling with SurPyval

Sometimes failure data is scarce or takes too long to collect: items are highly reliable, test time is limited, and few (or no) units fail during the observation window. But failure is often the end point of a gradual, measurable process — a crack grows, a resistance drifts, a lumen output fades, a material wears. Degradation analysis exploits this: instead of waiting for units to fail, we track a degradation measurement over time on each unit, define failure as the measurement crossing a threshold, and extrapolate each unit’s degradation trend to that threshold to obtain a failure time — even for units that never actually failed on test.

For the concepts and the underlying theory — the general-path model, the Lu-Meeker two-stage correction, the induced failure-time distribution, and the stochastic-process (Wiener/Gamma) first-passage models — see the Degradation Analysis page.

SurPyval implements the classic pseudo-failure-time approach:

1. A degradation path model (e.g. linear) is fitted, by least squares, to each unit’s measurements.

2. Each unit’s fitted path is extrapolated to the failure threshold. The crossing time is that unit’s pseudo failure time.

3. A lifetime distribution (Weibull by default) is fitted to the pseudo failure times, and can then be used like any other SurPyval parametric model.

If a unit’s fitted path never reaches the threshold (for example, the unit is not degrading, or is trending away from the threshold), the unit is treated as right censored at its last observed time, and the censoring is passed through to the lifetime distribution fit.

Degradation path models

The path models available, and the pseudo failure time each implies for a threshold \(y_{t}\), are:

Name	Path	Threshold crossing time
"linear"	\(y = a + b x\)	\((y_{t} - a) / b\)
"quadratic"	\(y = a + b x + c x^{2}\)	first positive root of \(c t^{2} + b t + (a - y_{t})\)
"exponential"	\(y = a e^{b x}\)	\(\ln(y_{t} / a) / b\)
"offset-exponential"	\(y = a + b e^{c x}\)	\(\ln((y_{t} - a) / b) / c\)
"power"	\(y = a x^{b}\)	\((y_{t} / a)^{1/b}\)
"logarithmic"	\(y = a + b \ln(x)\)	\(e^{(y_{t} - a) / b}\)
"lloyd-lipow"	\(y = a - b / x\)	\(b / (a - y_{t})\)
"gompertz"	\(y = a e^{-b e^{-c x}}\)	\(-\ln(-\ln(y_{t}/a)/b) / c\)
"michaelis-menten"	\(y = a x / (b + x)\)	\(b y_{t} / (a - y_{t})\)

Degradation can be increasing (crack length) or decreasing (luminous flux); the direction is captured by the sign of the fitted parameters and needs no configuration. Models that are linear in their parameters (linear, quadratic, logarithmic, Lloyd-Lipow) are fitted in closed form; the others are fitted by nonlinear least squares started from a linearised fit. The offset-exponential covers growth or decay toward an asymptote (a = 0 reduces it to the exponential); Gompertz is S-shaped; Michaelis-Menten saturates from zero toward a.

Not sure which shape fits? Pass path="best":

model = DegradationAnalysis.fit(x, y, i, threshold=150, path="best")
model.path_model.name   # the selected model
model.path_selection    # AICc score per candidate

Every registered path model is fitted to every unit and the model with the smallest AICc (pooled over all units, penalising the per-unit parameter count) is selected; candidates that cannot be fitted to every unit — domain violations such as negative measurements for the exponential, too few distinct measurement times for their parameter count, or non-convergence — are excluded and score nan. The selection is by measurement fit only; as always, prefer a shape with physical justification when one is known, since the winner is extrapolated well beyond the data.

Example

Consider four units whose degradation is measured every 100 hours, with failure defined as the measurement reaching 150:

import numpy as np
from surpyval.degradation import DegradationAnalysis

x = np.tile(np.arange(100, 1100, 100), 4)
i = np.repeat([1, 2, 3, 4], 10)
slopes = np.repeat([0.31, 0.28, 0.44, 0.37], 10)
y = 10 + slopes * x

model = DegradationAnalysis.fit(x, y, i, threshold=150)
print(model)

Degradation Analysis SurPyval Model
===================================
Path Model          : Linear
Threshold           : 150.0
Number of Units     : 4
Censored Units      : 0
Life Distribution   : Weibull
Parameters          :
     alpha: 441.4780882117898
      beta: 6.987078993008337

The fitted model exposes the per-unit results and forwards the usual lifetime functions to the fitted life model:

model.pseudo_failure_times
# array([451.61290323, 500.        , 318.18181818, 378.37837838])

model.sf([300, 400, 500])
# array([0.93496716, 0.60538471, 0.09196813])

model.life_model    # the underlying Parametric Weibull model
model.plot()        # data, fitted paths, and the threshold

Data can also come straight from a DataFrame with DegradationAnalysis.fit_from_df(df, x="time", y="measurement", i="unit", threshold=150), the life distribution and fitting method can be changed with the distribution and how arguments (e.g. distribution=LogNormal, how="MPP"), and a custom path shape can be used by passing a surpyval.degradation.PathModel subclass instance as path.

Predicting a new unit’s failure time

A fitted model can estimate the failure time of a new, partially observed unit from its degradation trajectory. The model fits its path shape to the new measurements and extrapolates to the same threshold:

# a new unit observed for only 300 hours, degrading at ~0.35/hour
x_new = [100, 200, 300]
y_new = [45, 80, 115]

model.predict_failure_time(x_new, y_new)   # ~400: when y reaches 150
model.predict_remaining_life(x_new, y_new) # ~100: minus its age (300)

If the trajectory has already crossed the threshold, the predicted failure time is in the past and the remaining life is negative. If the new unit’s fitted path never reaches the threshold (it is not degrading), both return nan with a warning. For a population-level view instead of a per-unit extrapolation, the fitted life model can be used directly — e.g. the survival of a unit that has already survived to time a: model.life_model.cs(t, a).

Bayesian remaining-life prediction

predict_failure_time trusts the new unit’s least-squares fit completely — dangerous when the trajectory is short or noisy. predict_rul instead blends the unit’s own trend with the population: the population path-parameter distribution (below) is the prior, the unit’s measurements are the likelihood, and the Gaussian posterior of the unit’s path parameters is pushed through the threshold crossing by Monte Carlo:

pred = model.predict_rul(x_new, y_new, alpha_ci=0.05)

pred.failure_time           # posterior median failure time
pred.failure_time_interval  # 95% credible interval
pred.rul                    # median remaining useful life
pred.rul_interval           # 95% credible interval
pred.prob_failed            # P(already crossed the threshold)
pred.prob_never_fails       # P(path never reaches the threshold)
pred.posterior_mean         # the unit's posterior path parameters
pred.posterior_cov

The posterior mean is a precision-weighted compromise: a short or noisy trajectory is shrunk toward the population’s typical path, and as measurements accumulate the prediction converges to the plain least-squares extrapolation. Unlike predict_failure_time, it works from a single measurement, and a not-yet-degrading trajectory yields a long-but-finite prediction with wide bounds rather than nan. The posterior is exact (conjugate) for path models that are linear in their parameters (linear, quadratic, logarithmic, Lloyd-Lipow) and an iterated-linearisation (Laplace) approximation for the others. It requires a positive measurement_var — with noiseless training paths there is nothing to blend.

The population path-parameter distribution

The fitted model also estimates the population distribution of the path parameters, which is what a random-effects treatment (and any Bayesian blending of a new unit’s trajectory with the population) needs as its prior:

model.path_param_mean        # mean path parameters across units
model.path_param_cov         # between-unit covariance (corrected)
model.path_param_sample_cov  # raw sample covariance (uncorrected)
model.measurement_var        # pooled measurement-error variance

Because each unit’s fitted parameters are least-squares estimates, their scatter across units mixes two sources: real unit-to-unit variability and per-unit estimation noise (\(\mathrm{Cov}(\hat{\theta}_i) = \Sigma + V_i\)). The raw sample covariance therefore overstates the between-unit variability. path_param_cov applies the Lu-Meeker two-stage correction: the measurement variance is pooled from the per-unit residuals, each unit’s estimation covariance \(V_i = \sigma^2 (J_i^T J_i)^{-1}\) is computed from the path Jacobian, and the average is subtracted from the sample covariance. The result is projected onto the positive semi-definite cone; if material clipping was needed (estimation noise comparable to the between-unit scatter — few units or few measurements per unit), a warning is raised and the corrected covariance should be treated as unreliable. When every unit has only as many measurements as path parameters, the measurement variance cannot be estimated and no correction is applied.

REML estimation of the population

The moments correction can go rank-deficient when the estimation noise rivals the between-unit scatter. The robust alternative is to fit the random-effects (Lu-Meeker) formulation directly as a linear mixed model — each unit’s parameters are draws \(\theta_i \sim MVN(\mu, \Sigma)\), so with the random effects integrated out each unit’s measurement vector is marginally

\[y_i \sim N(X_i \mu, \; X_i \Sigma X_i^T + \sigma^2 I)\]

and \((\mu, \Sigma, \sigma^2)\) are estimated by maximising the restricted (REML) marginal likelihood — REML rather than plain ML so the variance components do not inherit the small-sample downward bias from estimating \(\mu\). Select it with:

model = DegradationAnalysis.fit(
    x, y, i, threshold=150, population_method="reml"
)

The estimates land in the same attributes (path_param_mean, path_param_cov, measurement_var), so predict_rul and everything else work unchanged. \(\Sigma\) is parameterised by its Cholesky factor, so it is positive definite by construction — no clipping. On balanced designs (every unit measured at the same times) REML coincides with the corrected moments estimate; they differ on unbalanced data and when the unit count is small, where REML is preferable. REML requires a positive measurement variance.

For path models that are linear in their parameters (linear, quadratic, logarithmic, Lloyd-Lipow) the design matrix \(X_i\) is fixed and the marginal model above is exact. For a nonlinear path (exponential, power, Gompertz, …) the mean \(f(x_i, \theta_i)\) is no longer \(X_i \theta_i\), so REML uses the Lindstrom–Bates FOCE linearisation [LindstromBates1990]: each unit’s parameters are estimated at their conditional (penalised-least-squares) mode, the path is linearised about that mode to give a working linear mixed model, and the linear REML step is iterated to convergence. On a linear path this reduces to the exact fit in a single pass. Select it the same way:

model = DegradationAnalysis.fit(
    x, y, i, threshold=20, path="exponential",
    population_method="reml",
)

Induced failure-time distribution (Lu-Meeker)

Everything so far follows the pseudo-failure-time route: extrapolate each unit’s fitted path to the threshold, get one (noisy) failure time per unit, and fit a lifetime distribution to those times. That is simple and robust, but it throws away structure — it treats each extrapolated time as a plain observation and never uses the fact that we have estimated the whole population of paths.

The Lu-Meeker approach uses that population directly. Once the model has fitted the population path-parameter distribution \(\theta \sim N(\mu, \Sigma)\) (the path_param_mean and path_param_cov above), the failure-time distribution of the population follows deterministically: a unit with path parameters \(\theta\) fails at the time its path crosses the threshold, i.e. at \(T(\theta) = \text{inv\_path}(D; \theta)\). So the population failure-time distribution is simply the distribution of \(T(\theta)\) as \(\theta\) ranges over its population. There is rarely a closed form, so induced_life evaluates it by Monte Carlo: draw many \(\theta \sim N(\mu, \Sigma)\), push each through inv_path, and collect the resulting failure times.

import numpy as np
import matplotlib.pyplot as plt
from surpyval.degradation import DegradationAnalysis

rng = np.random.default_rng(0)
threshold = 30.0

xs, ys, ids = [], [], []
for unit in range(50):
    a = rng.normal(0.0, 0.2)      # per-unit intercept
    b = rng.normal(1.0, 0.25)     # per-unit slope (degradation rate)
    t = np.arange(0, 20, 2.0)
    y = a + b * t + rng.normal(0, 0.4, t.size)
    xs.append(t)
    ys.append(y)
    ids.append(np.full(t.size, unit))
x, y, i = (np.concatenate(z) for z in (xs, ys, ids))

model = DegradationAnalysis.fit(x, y, i, threshold=threshold, path="linear")
induced = model.induced_life(n_samples=20000, random_state=1)
induced

InducedFailureDistribution(Linear path, threshold=30, median=28.6674, prob_never_fails=0)

The returned InducedFailureDistribution behaves like any other life object — sf, ff, qf, mean, median and random are all there — so you can read life quantiles straight off it:

print("pseudo-failure median :", round(float(model.qf(0.5)), 2))
print("induced median        :", round(induced.median(), 2))
print("induced B10 life      :", round(induced.qf(0.10), 2))

pseudo-failure median : 30.29
induced median        : 28.67
induced B10 life      : 22.56

The real purpose is as a diagnostic. The pseudo-failure fit and the induced distribution come at the population life two different ways; when the path model and its population summary are trustworthy, they should agree. Overlaying the two CDFs is the check:

t = np.linspace(0, 60, 200)
plt.plot(t, model.ff(t), label="pseudo-failure fit")
plt.plot(t, induced.ff(t), "--", label="induced (Lu-Meeker)")
plt.xlabel("Time")
plt.ylabel("Probability of failure  F(t)")
plt.legend()

<matplotlib.legend.Legend at 0x7812f15cf6b0>

Close agreement (as here) is reassuring; a large gap is a warning that the path model, the Gaussian population assumption, or the covariance estimate is off.

A subtlety the induced distribution surfaces honestly: some draws of \(\theta\) describe paths that never reach the threshold (a non-increasing slope, say). Those contribute an inf failure time — a defective “never fails” mass reported as prob_never_fails — and once the quantiles reach into that mass they, and the mean, become inf. This is the correct behaviour: if a fraction of the population genuinely never fails, the population has no finite mean life. Finally, induced_life is defined for the plain (non-accelerated) population; an accelerated (covariate) model’s path parameters are not yet stress-conditional, so it is refused there.

Confidence bounds

The pseudo-failure-time approach is a two-stage estimator: the life distribution is fitted to extrapolated failure times as if they had been observed exactly. The plain life-model bounds therefore treat the pseudo failure times as certain and are too narrow. DegradationModel.cb corrects this, folding the first-stage (path-fit and extrapolation) uncertainty back into the life-model covariance with an analytic delta-method / generated-regressor correction:

import numpy as np
from matplotlib import pyplot as plt
from surpyval.degradation import DegradationAnalysis

rng = np.random.default_rng(0)
times = np.arange(1, 9) * 100.0
n_units = 60
x = np.tile(times, n_units)
unit = np.repeat(np.arange(n_units), times.size)
slopes = rng.normal(0.22, 0.05, size=n_units)     # between-unit spread
y = 10 + np.repeat(slopes, times.size) * x + rng.normal(0, 2.0, size=x.size)

model = DegradationAnalysis.fit(x, y, unit, threshold=150)

t = np.linspace(400, 800, 200)
band = model.cb(t, on='sf')                # (n, 2): two-stage [lower, upper]
plt.plot(t, model.sf(t), 'b', label='S(t)')
plt.fill_between(t, band[:, 0], band[:, 1], alpha=0.2,
                 label='95% two-stage band')
plt.legend()
plt.xlabel('Time')
plt.ylabel('S(t)')

Text(0, 0.5, 'S(t)')

The correction adds a positive term to the life-model information inverse, so the two-stage parameter covariance (model.life_parameter_covariance()) is the ordinary MLE covariance plus the propagated first-stage variance — the bounds widen to their correct coverage. A slower, assumption-light cross-check resamples whole units and reruns the whole pipeline:

model.cb(np.array([500.0, 600.0]), on='sf', method='bootstrap',
         n_boot=100, seed=0)

array([[0.739197  , 0.88378125],
       [0.54991711, 0.6958987 ]])

Random-effects (Lu-Meeker, fitted by REML) and stochastic-process (Wiener, gamma-process) degradation models, which propagate this uncertainty through a single likelihood rather than a two-stage correction, are candidates for future work.

Accelerated degradation testing (covariates)

In accelerated degradation testing (ADT) units are run at elevated stress (temperature, voltage, load) so they degrade fast enough to measure, and life is then extrapolated back to use conditions. Passing a per-unit stress covariate Z to DegradationAnalysis.fit() fits a regression life model on the pseudo failure times instead of a plain distribution, so life can be predicted at any stress. A plain distribution is wrapped automatically in an accelerated-failure-time model; an explicit regression fitter (AFT(LogNormal), WeibullPH, …) is used as given.

rng = np.random.default_rng(0)
times = np.arange(1, 11) * 5.0
xs, ys, ids, Zs = [], [], [], []
uid = 0
for stress in [0.0, 0.5, 1.0, 1.5]:        # four stress levels
    for _ in range(12):
        rate = 0.5 * np.exp(0.8 * stress) * np.exp(rng.normal(0, 0.1))
        path = 10 + rng.normal(0, 1) + rate * times
        xs.append(times)
        ys.append(path + rng.normal(0, 0.5, times.size))
        ids.append(np.full(times.size, uid))
        Zs.append(np.full(times.size, stress))
        uid += 1
xd, yd, idd, Zd = (np.concatenate(a) for a in (xs, ys, ids, Zs))

model = DegradationAnalysis.fit(xd, yd, idd, threshold=100.0, Z=Zd)
model

Degradation Analysis SurPyval Model
===================================
Path Model          : Linear
Threshold           : 100.0
Number of Units     : 48
Censored Units      : 0
Life Distribution   : Weibull (Log Linear [exp(beta'Z)] covariates)
Parameters          :
     alpha: 191.36340888317648
      beta: 10.908555487738695
    beta_0: 0.8513343841925708

The last fitted coefficient is the stress effect (higher stress ⇒ faster degradation ⇒ shorter life). The prediction methods now take the stress vector Z at which to evaluate life, so life at use conditions is one call:

for stress in [0.0, 0.5, 1.0]:
    print(f'stress {stress}: mean life = {model.mean(Z=[stress]):.1f}')

t = np.linspace(0, 300, 200)
for stress in [0.0, 0.5, 1.0]:
    plt.plot(t, model.sf(t, Z=[stress]), label=f'stress {stress}')
plt.legend()
plt.xlabel('Time')
plt.ylabel('Reliability at stress')

stress 0.0: mean life = 182.7
stress 0.5: mean life = 119.4
stress 1.0: mean life = 78.0

Text(0, 0.5, 'Reliability at stress')

qf and mean invert / integrate the regression survival function, and random draws from it.

Two-stage confidence bounds at a stress are available by bootstrap: units are resampled (each carrying its stress), the whole accelerated pipeline is rerun, and the reliability at Z is read off each refit, so the first-stage path/extrapolation uncertainty is folded in — just as for the plain model, but evaluated at a chosen stress:

t = np.array([50.0, 100.0, 150.0])
band = model.cb(t, on='sf', method='bootstrap', Z=[0.0],
                n_boot=100, seed=0)
band                                        # (n, 2): [lower, upper] at Z=0

array([[0.99999697, 0.99999999],
       [0.99754806, 0.9998559 ],
       [0.86674985, 0.97136315]])

The analytic (generated-regressor) delta-method correction used for the plain model is not derived for the regression life fit, so method='bootstrap' is required for a covariate model (and model.cb needs the stress Z). The first-stage-only regression bounds — which ignore the extrapolation uncertainty — remain available directly through model.life_model.cb(x, Z, ...).

Stochastic-process degradation models

Everything above is the general-path approach: fit a deterministic curve to each unit, extrapolate it to the threshold to get a pseudo failure time, then fit a lifetime distribution to those times. It works well when each unit really does follow a smooth trend plus measurement noise.

But sometimes the degradation itself is random over time, not a smooth curve observed with error. A crack does not grow along a tidy exponential; it jumps ahead in fits and starts. A wear signal wanders up and down from measurement to measurement. In those cases it is more honest to model the increments of the degradation as a stochastic process, and read the failure-time distribution off the process directly — as the distribution of the first time the process crosses the threshold (the “first-passage time”). This is a standard part of the reliability toolkit; see [Meeker1998] for a textbook treatment and [LuMeeker1993] for the idea of deriving a failure-time distribution from degradation measurements.

SurPyval provides two such processes. They are not competitors; they describe different physics, and the right one is dictated by whether your degradation can decrease:

• WienerProcess — for signals that fluctuate up and down (noisy sensors, measurements that wobble).

• GammaProcess — for damage that only ever accumulates (wear, corrosion, crack growth).

Both are fitted from the same three arrays you have used throughout this section: x (measurement times), y (degradation measurements) and i (the unit each measurement belongs to). Internally each model looks only at the increments between consecutive measurements of a unit, so units can be measured at different, irregular times without any special handling — a two-week gap simply contributes a larger dt.

The Wiener process

A Wiener process with drift (Brownian motion with drift) models the degradation as

\[W(t) = \mu\, t + \sigma\, B(t),\]

where \(B(t)\) is standard Brownian motion. There are just two parameters, and it is worth being clear about what each one means:

• \(\mu\) — the drift. This is the average rate at which degradation accumulates: on average the signal climbs by \(\mu\) per unit time. A larger drift means faster wear-out and a shorter life.

• \(\sigma\) — the diffusion (or volatility). This is the size of the random wobble around that average trend. With \(\sigma = 0\) the process would be a perfectly straight line \(\mu t\); the bigger \(\sigma\), the more the path jitters up and down and the more spread-out the failure times become. (Wiener degradation models, including extensions with unit-to-unit random effects, are surveyed in [Wang2010].)

Over any interval of length \(\Delta t\), the change in degradation is Gaussian:

\[\Delta W \sim \mathrm{Normal}\!\left(\mu\, \Delta t,\; \sigma^2\, \Delta t\right).\]

Two things follow from this. First, because a Normal can be negative, the path can go down as well as up — which is exactly why the Wiener process is the right model for noisy, non-monotone signals. Second, the increments are independent, so fitting is easy: the drift is just the total degradation divided by the total time, and the diffusion is estimated from how much the increments scatter around that trend. SurPyval does this by maximum likelihood.

Failure = first passage. A unit fails the first time \(W(t)\) reaches the threshold \(D\). For a Wiener process this first-passage time has a famous closed form — the Inverse Gaussian distribution — with

\[\text{mean} = \frac{D}{\mu}, \qquad \text{shape} = \frac{D^2}{\sigma^2}.\]

The mean life \(D/\mu\) is beautifully intuitive: distance to failure divided by the average speed. The shape controls how tightly the failure times cluster around that mean (more diffusion → more scatter). Because the increments are independent, the model also needs the drift to be positive — a non-positive drift would mean the process is not reliably heading toward the threshold at all, so SurPyval raises an error rather than return a “life” that may never end.

Let’s fit one. We simulate 30 units, each measured every half–time-unit, with a true drift of 0.5 and diffusion 0.4, failing at a degradation of 10:

import numpy as np
import matplotlib.pyplot as plt
from surpyval.degradation import WienerProcess

rng = np.random.default_rng(0)
mu_true, sigma_true, threshold = 0.5, 0.4, 10.0

xs, ys, ids = [], [], []
for unit in range(30):
    t = np.arange(0, 15.5, 0.5)
    increments = rng.normal(
        mu_true * 0.5, sigma_true * np.sqrt(0.5), size=t.size - 1
    )
    y = np.concatenate([[0.0], np.cumsum(increments)])
    xs.append(t)
    ys.append(y)
    ids.append(np.full(t.size, unit))
x, y, i = (np.concatenate(a) for a in (xs, ys, ids))

model = WienerProcess.fit(x, y, i, threshold=threshold)
model

Wiener Process Degradation Model
================================
Drift (mu)          : 0.486346
Diffusion (sigma)   : 0.389563
Threshold           : 10
Mean time to failure: 20.5615

Read the summary line by line: the fitted drift and diffusion are close to the 0.5 and 0.4 we simulated, and the mean time to failure is threshold / drift — roughly 10 / 0.5 = 20 time units. Notice the paths below are jagged and occasionally dip downward; that non-monotone wobble is the Wiener process’s defining feature.

for unit in range(8):
    m = i == unit
    plt.plot(x[m], y[m], alpha=0.6)
plt.axhline(threshold, color="k", linestyle="--", label="threshold")
plt.xlabel("Time")
plt.ylabel("Degradation")
plt.legend()

<matplotlib.legend.Legend at 0x7812f0c25af0>

Now the payoff: a full failure-time distribution, derived from the process, that you can query like any other SurPyval model. mean() is the average life, ff(t) is the probability of having failed by time t (the CDF), sf(t) is the reliability, and qf(p) is the quantile (e.g. the median life, or the time by which 10 % have failed):

print("mean life           :", round(model.mean(), 2))
print("P(fail by t = 25)   :", round(model.ff(25.0), 3))
print("median life         :", round(model.qf(0.5), 2))
print("B10 life (10% fail) :", round(model.qf(0.10), 2))

t = np.linspace(0, 40, 200)
plt.plot(t, model.ff(t))
plt.xlabel("Time to failure")
plt.ylabel("Probability of failure  F(t)")

mean life           : 20.56
P(fail by t = 25)   : 0.885
median life         : 20.25
B10 life (10% fail) : 16.17

Text(0, 0.5, 'Probability of failure  F(t)')

Remaining useful life. The real power of a degradation model is that it can update its forecast for a unit you have been watching. If a unit is currently at degradation level 7 (out of a threshold of 10), only the remaining distance of 3 matters, and — because Wiener increments are independent of the past — the remaining life is itself an Inverse Gaussian over that shorter distance. predict_rul returns its median and an interval:

rul = model.predict_rul(current_degradation=7.0)
print("median remaining life :", round(rul.rul, 2))
print("90% interval          :", tuple(round(v, 2) for v in rul.rul_interval))

median remaining life : 5.87
90% interval          : (3.18, 10.88)

The interval widths tell you how much uncertainty remains: a unit close to the threshold has a short, tight remaining-life estimate; a fresh unit has a long, uncertain one. If the current degradation is already at or beyond the threshold, prob_already_failed is 1 and the remaining life is 0.

The Gamma process

The Wiener process allows the signal to decrease, which is wrong for damage that is physically irreversible — a crack never heals, corrosion never reverses, wear never un-wears. For those, use a Gamma process (whose use in maintenance and reliability is surveyed in [vanNoortwijk2009]), whose increments are strictly non-negative, so the path is monotone increasing.

Over an interval of length \(\Delta t\), the Gamma-process increment is Gamma-distributed:

\[\Delta W \sim \mathrm{Gamma}\!\left(\text{shape} = \alpha\, \Delta t,\; \text{rate} = \beta\right).\]

Again, two parameters, and again it is worth knowing what they mean:

• \(\alpha\) — the shape rate. It controls how quickly and how steadily damage accrues. The amount of shape accumulated by time \(t\) is \(\alpha t\); a large \(\alpha\) gives many small, regular increments (a smooth-looking climb), a small \(\alpha\) gives fewer, larger, jumpier increments.

• \(\beta\) — the rate parameter of those increments. It sets their scale: the mean degradation accumulated per unit time is \(\alpha / \beta\), and the variance per unit time is \(\alpha / \beta^2\).

So \(\alpha/\beta\) is the Gamma process’s analogue of the Wiener drift — the average degradation speed — while \(\alpha\) alone governs how regular versus jumpy the accumulation is.

Failure = first passage, again, but now the monotonicity makes it especially clean: the process has crossed the threshold \(D\) by time \(t\) exactly when its level \(W(t)\) is at or above \(D\). So the probability of having failed by \(t\) is

\[F(t) = \Pr\!\left(W(t) \ge D\right),\]

which SurPyval evaluates with the (regularised) incomplete gamma function. There is no simpler closed form than that, but every method you need — sf, ff, df, qf, mean, random — is computed from it.

Here we simulate 40 units of monotone wear with shape rate 3 and rate 1.5 (so mean degradation speed \(\alpha/\beta = 2\) per unit time), failing at 30:

from surpyval.degradation import GammaProcess

rng = np.random.default_rng(1)
alpha_true, beta_true, threshold = 3.0, 1.5, 30.0

xs, ys, ids = [], [], []
for unit in range(40):
    t = np.arange(0, 12.5, 0.5)
    increments = rng.gamma(alpha_true * 0.5, 1.0 / beta_true, size=t.size - 1)
    y = np.concatenate([[0.0], np.cumsum(increments)])
    xs.append(t)
    ys.append(y)
    ids.append(np.full(t.size, unit))
x, y, i = (np.concatenate(a) for a in (xs, ys, ids))

model = GammaProcess.fit(x, y, i, threshold=threshold)
model

Gamma Process Degradation Model
===============================
Shape rate (alpha)  : 2.97309
Rate (beta)         : 1.50409
Threshold           : 30
Mean time to failure: 15.3452

The fitted alpha and beta recover the 3 and 1.5 we used, and the mean time to failure is about threshold / (alpha/beta) = 30 / 2 = 15. The paths this time only ever climb — no downward wobble is possible:

for unit in range(8):
    m = i == unit
    plt.plot(x[m], y[m], alpha=0.6)
plt.axhline(threshold, color="k", linestyle="--", label="threshold")
plt.xlabel("Time")
plt.ylabel("Degradation")
plt.legend()

<matplotlib.legend.Legend at 0x7812ee968440>

The failure-time distribution and remaining-life prediction work exactly as they did for the Wiener model — same method names, same meaning:

print("mean life      :", round(model.mean(), 2))
print("median life    :", round(model.qf(0.5), 2))

rul = model.predict_rul(current_degradation=20.0)
print("RUL at y = 20  :", round(rul.rul, 2),
      "  interval", tuple(round(v, 2) for v in rul.rul_interval))

mean life      : 15.35
median life    : 15.29
RUL at y = 20  : 5.17   interval (2.85, 7.93)

Because a Gamma process cannot go down, passing it degradation that decreases over an interval is a modelling error, and SurPyval says so rather than fitting something meaningless — the message points you at the Wiener process instead:

# this data dips from 5 back to 3 -- not allowed for a monotone process
GammaProcess.fit([0, 1, 2], [0.0, 5.0, 3.0], [1, 1, 1], threshold=10.0)

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[16], line 2
      1 # this data dips from 5 back to 3 -- not allowed for a monotone process
----> 2 GammaProcess.fit([0, 1, 2], [0.0, 5.0, 3.0], [1, 1, 1], threshold=10.0)

File ~/checkouts/readthedocs.org/user_builds/surpyval/envs/latest/lib/python3.12/site-packages/surpyval/degradation/process_models.py:558, in GammaProcess.fit(cls, x, y, i, threshold)
    556 dt, dy = _increments(x, y, i)
    557 if np.any(dy < 0):
--> 558     raise ValueError(
    559         "the degradation decreases over at least one interval, but a "
    560         "Gamma process is monotone increasing. Use WienerProcess for "
    561         "non-monotone / noisy signals."
    562     )
    563 # Guard against exact-zero increments (log 0) by nudging.
    564 dy = np.where(dy <= 0, 1e-12, dy)

ValueError: the degradation decreases over at least one interval, but a Gamma process is monotone increasing. Use WienerProcess for non-monotone / noisy signals.

Choosing between Wiener and Gamma

The decision is almost always settled by one question: can your degradation physically decrease?

Your degradation paths…	Use	Because
fluctuate up and down (noisy sensor, wandering signal)	WienerProcess	Gaussian increments allow decreases and absorb measurement noise; failure time is a closed-form Inverse Gaussian.
only ever increase (wear, corrosion, crack growth, fatigue)	GammaProcess	non-negative increments enforce monotone, irreversible damage; forcing a Wiener fit would misread real jumps as noise.

A quick practical test: plot your paths. If they visibly dip, the Wiener process models that directly. If they are monotone, the Gamma process is the physically honest choice — and if you try to give monotone-only data that happens to dip (usually a measurement glitch) to the Gamma process, the error it raises is a useful signal in itself.

Compared with the general-path approach at the top of this page, both process models share two advantages: they handle irregular measurement spacing without any special treatment (each increment simply carries its own dt), and they give the failure-time distribution directly from the fitted process rather than through noisy per-unit pseudo failure times. The general-path models remain the better choice when each unit truly follows a smooth deterministic trend observed with error, or when you need a specific parametric path shape.

Saving and loading a fitted model

Every fitted degradation model can be serialised to a dictionary or JSON file and rebuilt later — the general-path DegradationModel, the stochastic- process WienerProcessModel / GammaProcessModel, and the InducedFailureDistribution:

from surpyval.degradation import DegradationAnalysis, DegradationModel

saveable = DegradationAnalysis.fit(
    np.tile(np.arange(100, 1100, 100), 4),
    10 + np.repeat([0.31, 0.28, 0.44, 0.37], 10)
    * np.tile(np.arange(100, 1100, 100), 4),
    np.repeat([1, 2, 3, 4], 10),
    threshold=150,
)
reloaded = DegradationModel.from_dict(saveable.to_dict())
grid = np.array([300.0, 450.0, 600.0])
print("match:", np.allclose(saveable.sf(grid), reloaded.sf(grid)))

match: True

Use to_json / from_json for a file directly. DegradationModel stores its raw data, so the reloaded model can also produce bootstrap confidence bounds; its fitted life model (plain or accelerated) round-trips through its own serialisation.