Degradation Analysis

Models for measurements of a unit’s condition over time – wear, crack length, capacity loss – from which the time to failure is inferred as the time the degradation crosses a threshold. Import everything on this page from surpyval.degradation. There are three approaches: a path model fitted to each unit’s measurements (DegradationAnalysis, which returns a DegradationModel); a stochastic process (the Wiener and gamma processes); and destructive degradation, where each unit is measured once (DestructiveDegradation). The theory is in Degradation Analysis and worked examples are in Degradation Modelling with SurPyval.

Degradation Analysis Fitter

DegradationAnalysis is an instance of the class below.

class surpyval.degradation.degradation_analysis.DegradationAnalysis_

Bases: object

Pseudo-failure-time degradation analysis.

Fits a degradation path model to each unit’s measurements, extrapolates each fitted path to the failure threshold to obtain per-unit pseudo failure times, and fits a lifetime distribution to those times. Units whose fitted path never reaches the threshold at a positive finite time are right censored at their last observed time (with a warning); units already past the threshold at their first measurement – the path crossed at or before time zero – have failed by then and are left censored at their first measurement time (with a warning). The side of the threshold that counts as failed is read from the units whose paths cross it.

Examples

>>> import numpy as np
>>> from surpyval.degradation import DegradationAnalysis
>>> # 4 units measured every 100 hours; degradation grows linearly
>>> # at a different rate per unit; failure is defined at level 150.
>>> x = np.tile(np.arange(100, 1100, 100), 4)
>>> slopes = np.repeat([0.31, 0.28, 0.44, 0.37], 10)
>>> i = np.repeat([1, 2, 3, 4], 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.47809611105606
      beta: 6.987078889297555
>>> model.pseudo_failure_times
array([451.61290323, 500.        , 318.18181818, 378.37837838])

fit(x: ArrayLike, y: ArrayLike, i: ArrayLike, threshold: float, path: str | ~surpyval.degradation.path_models.PathModel = 'linear', distribution: ~typing.Any = <surpyval.univariate.parametric.distributions.weibull.Weibull_ object>, how: str = 'MLE', population_method: str = 'moments', Z: ArrayLike | None = None, links: dict[str, str] | None = None, acceleration: str | None = None, stress_ref: ArrayLike | None = None) → DegradationModel

Fit a degradation analysis model.

Parameters:

x (array like) – Times at which the degradation measurements were taken.
y (array like) – The degradation measurements.
i (array like) – The unit each measurement belongs to. Must have the same length as x and y.
threshold (float) – The degradation level at which a unit is defined to have failed.
path (str or PathModel, optional) – The degradation path model fitted to each unit: one of "linear" (default), "quadratic", "exponential", "offset-exponential", "power", "logarithmic", "lloyd-lipow", "gompertz", "michaelis-menten", a PathModel instance, or "best" to fit every registered model to all units and select the one with the smallest AICc (the per-candidate scores are exposed as path_selection on the returned model; candidates that cannot be fitted to every unit are excluded).
distribution (ParametricFitter, optional) – The lifetime distribution fitted to the pseudo failure times. Defaults to Weibull.
how (str, optional) – The method used to fit the lifetime distribution (passed to distribution.fit). Defaults to "MLE".
population_method (str, optional) – How the population path-parameter distribution (path_param_mean, path_param_cov, measurement_var) is estimated. "moments" (default) uses the two-stage noise-corrected sample moments; "reml" maximises the restricted marginal likelihood of the mixed model, which never needs clipping and is preferable with few units; both warn when the estimated covariance is singular (on the boundary: a variance of zero or a correlation of +-1). Linear-in-parameter paths (linear, quadratic, logarithmic, lloyd-lipow) are fitted as an exact linear mixed model; nonlinear paths (exponential, power, gompertz, …) are fitted by the Lindstrom-Bates FOCE linearisation. Either way REML requires measurement noise (some unit with more measurements than path parameters).
Z (array like, optional) – Stress covariates for accelerated degradation testing (ADT). When given, the life model is fitted as a regression on the pseudo failure times – log(pseudo failure time) = f(Z) + noise – so that life can be predicted at any stress. Z is aligned to x/y/i (one row per measurement) and must be constant within each unit (a unit is tested at a single stress); it is reduced to one covariate row per unit. If distribution is already a regression fitter (e.g. AFT(Weibull), WeibullPH, CoxPH) it is used directly; a plain distribution (e.g. Weibull) is wrapped in an accelerated failure time model, AFT(distribution). The returned model’s prediction methods (sf, ff, qf, random …) then take the stress vector Z at which to evaluate life.
links (dict, optional) – Model the degradation mechanism against stress as well (requires Z): the path parameters named here depend on the unit’s stress, the rest do not. Each value is the link the parameter is modelled on – "identity" (the parameter is linear in Z) or "log" (its log is linear in Z, so a rate with Z = 1/T follows an Arrhenius relationship and the parameter stays positive). Per unit, eta_i = D(z_i) gamma + u_i on the link scale with a between-unit random effect u_i ~ MVN(0, Sigma); gamma and Sigma are estimated by the same two-stage or REML route as the plain population, and stored as path_param_fixed (labelled by path_param_fixed_names) and path_param_link_cov. The life model is still the covariate regression on the pseudo failure times, so every prediction method works as without links. For example links={"b": "log"} with the linear path lets the degradation rate b accelerate log-linearly with stress while the intercept a (the initial state) is common.
acceleration ({None, "clock"}, optional) – "clock" models stress as speeding up the clock of every unit’s path, which allows Z to change during a unit’s test (a step-stress test) as well as between units. A unit at stress z ages AF(z) = exp(gamma' (z - stress_ref)) times faster than at the reference stress, and its path is the path model evaluated on the reference-stress time it has aged, tau(t) = integral of AF(z(s)) ds. Z is then one row per measurement giving the stress applied over the interval that ends at that measurement (the first interval starts at time zero, so times must be non-negative). The path parameters, their population and the pseudo failure times are all on the reference-stress clock; distribution is fitted to those reference-stress lifetimes, and the prediction methods take the stress as Z (one stress row or a StepSchedule) to give life under any stress history, F(t) = F0(tau(t)). The stress coefficients are stored as gamma. With population_method="moments" they are estimated by profile least squares, which needs units whose stress changes during the test (a unit held at one stress can absorb any acceleration into its own path parameters); with "reml" by the mixed model, which also uses the differences between units at different stresses. Cannot be combined with links or path="best".
stress_ref (array like, optional) – The reference stress for acceleration="clock" (usually the use condition), one row. Defaults to the mean stress over the measurement intervals.

Returns:

The fitted degradation model, with the per-unit paths, pseudo failure times, and the fitted life model.

Return type:

DegradationModel

fit_from_df(df: DataFrame, x_col: str = 'x', y_col: str = 'y', i_col: str = 'i', Z_cols: str | list[str] | None = None, **fit_kwargs: Any) → DegradationModel

Fit a degradation analysis model from a DataFrame.

The column arguments end in _col (_cols for a list), as in every fit_from_df (principle 21); their v0.21 names x, y and i still work, with a DeprecationWarning, until v0.23.

Parameters:

df (DataFrame) – DataFrame with the degradation data.
x_col (str, optional) – Column of the measurement times. Defaults to "x".
y_col (str, optional) – Column of the degradation measurements. Defaults to "y".
i_col (str, optional) – Column of the unit identifiers. Defaults to "i".
Z_cols (str or list of str, optional) – Column(s) of the stress covariates for accelerated degradation testing. When given, the selected columns are passed as Z to fit(), fitting a covariate (ADT) life model. Their names are recorded on the model as Z_cols (and kept by to_dict), so every method that takes Z also takes a DataFrame and selects these columns by name.
**fit_kwargs – Remaining arguments passed to fit(): threshold (required), and optionally path, distribution, how, population_method, links, acceleration and stress_ref.

Returns:

The fitted degradation model.

Return type:

DegradationModel

Degradation Model

class surpyval.degradation.degradation_analysis.DegradationModel(x: NDArray, y: NDArray, i: NDArray, units: NDArray, threshold: float, path_model: Any, path_params: NDArray, pseudo_failure_times: NDArray, c: NDArray, life_model: Any, measurement_var: float, path_param_mean: NDArray, path_param_cov: NDArray, path_param_sample_cov: NDArray, population_method: str, path_selection: dict | None = None, Z: NDArray | None = None, links: dict[str, str] | None = None, path_param_fixed: NDArray | None = None, path_param_fixed_names: list[str] | None = None, path_param_link_cov: NDArray | None = None, acceleration: str | None = None, gamma: ArrayLike | None = None, stress_ref: ArrayLike | None = None, Z_cols: list[str] | None = None)

Bases: SerialisableMixin

A fitted degradation analysis model.

This is the model object returned by DegradationAnalysis.fit(). It holds the per-unit fitted degradation paths, the pseudo failure times extrapolated from them, and the lifetime distribution fitted to those pseudo failure times. The usual lifetime functions (sf, ff, df, hf, Hf, qf, mean, random) are forwarded to the fitted life model, and the failure time of a new, partially observed unit can be estimated from its trajectory with predict_failure_time() / predict_remaining_life().

Parameters:

x (ndarray) – The degradation data: measurement times, measurements, and the unit each measurement belongs to.
y (ndarray) – The degradation data: measurement times, measurements, and the unit each measurement belongs to.
i (ndarray) – The degradation data: measurement times, measurements, and the unit each measurement belongs to.
units (ndarray) – The distinct unit identifiers.
threshold (float) – The degradation level at which a unit is considered failed.
path_model (PathModel) – The degradation path model fitted to each unit.
path_params (ndarray) – Per-unit fitted path parameters, one row per entry of units.
pseudo_failure_times (ndarray) – Per-unit pseudo failure time: the extrapolated threshold crossing time; the unit’s last observed time for a unit whose path never reaches the threshold; its first (positive) measurement time for a unit already past the threshold there.
c (ndarray) – Per-unit censor flags: 0 where the fitted path crosses the threshold at a positive time, 1 (right censored) where it never reaches it, -1 (left censored: failed before its first measurement) where the path is already past the threshold at the first measurement, having crossed at or before time zero.
life_model (Parametric) – The lifetime distribution fitted to the pseudo failure times.
measurement_var (float) – Pooled estimate of the measurement-error variance around the per-unit paths (the per-unit residual sums of squares over the total residual degrees of freedom). Zero when every unit has exactly as many measurements as path parameters.
path_param_mean (ndarray) – Mean of the per-unit fitted path parameters: the estimated population mean path.
path_param_cov (ndarray) – Noise-corrected estimate of the between-unit covariance of the true path parameters (Lu-Meeker two-stage): the sample covariance of the per-unit estimates minus the average least-squares estimation covariance, projected onto the positive semi-definite cone.
path_param_sample_cov (ndarray) – The raw (uncorrected) sample covariance of the per-unit fitted path parameters. This overstates the between-unit variability because each per-unit estimate also carries least-squares estimation noise.
population_method (str) – How the population estimates (measurement_var, path_param_mean, path_param_cov) were obtained: "moments" (two-stage correction) or "reml".
path_selection (dict or None) – When fitted with path="best", the AICc score of every candidate path model (nan for candidates that could not be fitted to every unit); None otherwise. The fitted path_model is the candidate with the smallest score. The keys are the models’ display names ("Offset Exponential"), not the path= strings.
Z (ndarray or None) – The stresses of an accelerated model: one row per unit (aligned to units) for a model fitted with Z alone or with links; for a step-stress (acceleration="clock") model the stress rows as given, one per measurement (aligned to x). None for a model fitted without stress.
links (dict or None) – When the path parameters were modelled against stress (links given to DegradationAnalysis.fit()), the stress-dependent parameters and their links; None otherwise. With links the population of path parameters is stress-conditional, on the link scale: eta_i = D(z_i) gamma + u_i with u_i ~ MVN(0, Sigma) and theta_i = h(eta_i).
path_param_fixed (ndarray or None) – The fixed effects gamma of the stress-conditional population model, labelled by path_param_fixed_names: for every path parameter its link-scale intercept, followed (for the stress-dependent ones) by its coefficient on each covariate.
path_param_fixed_names (list of str or None) – Labels for path_param_fixed: the link-scale parameter name ("log(b)" for a log link) and "<name>:Z<j>" for the coefficient on covariate j.
path_param_link_cov (ndarray or None) – The between-unit covariance Sigma of the link-scale path parameters given the stress – the scatter left after the stress effect is removed, unlike the pooled path_param_cov which mixes the stress levels.
acceleration (str or None) – "clock" when stress was modelled as speeding up the clock of every unit’s path (acceleration="clock" in DegradationAnalysis.fit()); None otherwise. The path parameters, their population, the pseudo failure times and the life model are then all on the reference-stress clock, and Z holds the stress rows aligned to x.
gamma (ndarray or None) – The stress coefficients of the clock: a unit at stress z ages exp(gamma' (z - stress_ref)) times faster than at the reference stress.
stress_ref (ndarray or None) – The reference stress of the clock.
Z_cols (list of str or None) – The covariate columns of a model fitted with DegradationAnalysis.fit_from_df(): every method that takes Z then also takes a DataFrame and selects these columns by name. None for a model fitted from arrays, which refuses a DataFrame.

Examples

Eight units, each degrading linearly at its own rate and measured ten times, fail when the measurement reaches 450:

>>> import numpy as np
>>> from surpyval.degradation import DegradationAnalysis
>>> rng = np.random.default_rng(1)
>>> x = np.tile(np.arange(100.0, 1100.0, 100.0), 8)
>>> i = np.repeat(np.arange(8), 10)
>>> a = np.repeat(rng.normal(10.0, 3.0, 8), 10)
>>> b = np.repeat(rng.normal(0.3, 0.05, 8), 10)
>>> y = a + b * x + rng.normal(0, 3.0, x.size)
>>> model = DegradationAnalysis.fit(x, y, i, threshold=450)
>>> model.pseudo_failure_times.round(1)
array([1393.1, 1377.5, 1455.4, 1358.9, 1673.9, 1495.6, 1592.4, 1334. ])

A Weibull is fitted to those, and the lifetime functions use it:

>>> model.life_model.params.round(3)
array([1515.035,   12.785])
>>> model.sf([1200, 1500]).round(4)
array([0.9505, 0.4147])

Hf(x: ArrayLike, Z: Any = None) → NDArray: Cumulative hazard of the life model (Z for accelerated).

Z: NDArray | None: Per-unit covariates when fitted as an accelerated-degradation model (Z given to DegradationAnalysis.fit()); None otherwise.

acceleration_factor(Z: Any) → float

How much faster a unit ages at stress Z than at the reference stress: exp(gamma' (z - stress_ref)), for a model fitted with acceleration="clock".

A unit held at Z degrades along its path this many times faster, and a life at the reference stress divides by it to give the life at Z.

Parameters:: Z (array like or DataFrame) – One stress row (nan for a missing value gives nan); a one-row DataFrame for a model fitted with fit_from_df.

cb(x: ArrayLike, Z: Any = None, on: str = 'sf', alpha_ci: float = 0.05, bound: str = 'two-sided', method: str = 'analytic', n_boot: int = 200, random_state: int | None = None) → NDArray

Confidence bounds on the reliability of the fitted life model that account for the degradation analysis being a two-stage estimator.

The pseudo failure times are extrapolated per-unit path fits, not observed failures, so life_model.cb – which treats them as exact – gives intervals that are too narrow. These bounds fold the first-stage (measurement + extrapolation) uncertainty back in.

For an accelerated-degradation (covariate) model the bound is evaluated at a stress vector Z and only method='bootstrap' is available: the generated-regressor delta-method correction is not derived for the regression life fit, so the first-stage uncertainty is folded in by resampling units (each carrying its stress) and rerunning the whole accelerated pipeline.

Parameters:

x (array like) – Times at which to evaluate the bound(s).
Z (array like, optional) – Stress vector at which to evaluate the bound; required for an accelerated model, rejected for a plain one. For a step-stress (acceleration="clock") model it is one stress row or a StepSchedule, and only method='bootstrap' is available: units are resampled with their stress histories and the clock is re-estimated on each resample (with the model’s population_method, so a "reml" model’s bootstrap takes correspondingly longer).
on ({'sf', 'ff', 'Hf'}, optional) – The function to bound ('R' and 'F' are accepted as aliases of 'sf' and 'ff'). Default 'sf'.
alpha_ci (float, optional) – Total tail probability of the bound(s): a two-sided band has alpha_ci / 2 in each tail. Default 0.05 (a 95% band).
bound ({'two-sided', 'lower', 'upper'}, optional) – Two-sided bounds put [lower, upper] on the last axis.
method ({'analytic', 'bootstrap'}, optional) – 'analytic' (default) is a fast delta-method correction; 'bootstrap' resamples units and reruns the whole pipeline (a slower, assumption-light cross-check). Accelerated (covariate) models support 'bootstrap' only.
n_boot (int, optional) – Bootstrap resamples (method='bootstrap' only). Default 200.
random_state (int or numpy.random.Generator, optional) – Seed or generator for the bootstrap resampling. None (the default) seeds from numpy’s global RNG, so np.random.seed controls it.

Returns:

The confidence bound(s) on on at each x.

Return type:

numpy array

df(x: ArrayLike, Z: Any = None) → NDArray: Density of the fitted life model (Z for accelerated models).

ff(x: ArrayLike, Z: Any = None) → NDArray: CDF of the fitted life model (pass Z for accelerated models).

classmethod from_dict(model_dict: dict) → DegradationModel

Rebuild a degradation model from a to_dict() dictionary.

The path model is resolved by name (its PATH_MODELS key, or the display name that older dictionaries stored) and the life model by its own from_dict; both are restricted to the known types.

See also

to_dict, to_json, from_json

classmethod from_json(fp: str | PathLike) → Any: Load a model from a JSON file written by to_json(), or from the JSON text it returned (a string starting with {).

hf(x: ArrayLike, Z: Any = None) → NDArray: Hazard rate of the fitted life model (Z for accelerated).

induced_life(n_samples: int = 10000, *, Z: Any = None, random_state: int | None = None) → InducedFailureDistribution

The population failure-time distribution induced by the path model (the Lu-Meeker approach), as a Monte-Carlo diagnostic complement to the pseudo-failure-time life_model.

Path parameters are drawn from the fitted population distribution theta ~ N(path_param_mean, path_param_cov) and each draw is pushed through the path model’s inv_path(threshold) to a failure time. This derives the population life directly from the path model, rather than via each unit’s noisy extrapolated failure time. Overlaying the returned distribution’s ff on this model’s own ff is a check that the two agree.

Parameters:

n_samples (int, optional) – Number of Monte-Carlo path-parameter draws. Default 10000.
Z (array like, optional) – The stress to induce the life at. Required for a model whose path parameters were modelled against stress (fitted with links): the draws are then eta ~ N(D(z) gamma, Sigma) on the link scale, mapped through the links to path parameters. Refused for a model without links, unless it is a step-stress (acceleration="clock") model: then Z is required, as one stress row or a StepSchedule, and each draw’s reference-stress failure time is read along that stress’s clock. The returned distribution records a constant stress row as its stress; under a profile it records none.
random_state (int or numpy.random.Generator, optional) – Seed or generator for a reproducible result. None (the default) seeds from numpy’s global RNG, so np.random.seed controls it.

Returns:

The Monte-Carlo induced failure-time distribution.

Return type:

InducedFailureDistribution

Examples

The induced median next to the pseudo-failure fit’s:

>>> import numpy as np
>>> from surpyval.degradation import DegradationAnalysis
>>> rng = np.random.default_rng(1)
>>> x = np.tile(np.arange(100.0, 1100.0, 100.0), 8)
>>> i = np.repeat(np.arange(8), 10)
>>> a = np.repeat(rng.normal(10.0, 3.0, 8), 10)
>>> b = np.repeat(rng.normal(0.3, 0.05, 8), 10)
>>> y = a + b * x + rng.normal(0, 3.0, x.size)
>>> model = DegradationAnalysis.fit(x, y, i, threshold=450)
>>> induced = model.induced_life(random_state=0)
>>> round(induced.median()), round(float(model.qf(0.5)))
(1452, 1472)
>>> induced.prob_never_fails
0.0

property is_accelerated: bool

the life model is a covariate (ADT) regression model, or stress accelerates the clock (acceleration="clock").

Type:: True when life depends on stress

life_model: Any: A plain Parametric life model, or the regression model for an accelerated (covariate) fit.

life_parameter_covariance(method: str = 'analytic') → NDArray

Covariance of the fitted life-model parameters, corrected for the first-stage (path-fit and extrapolation) uncertainty that the plain life-model MLE ignores.

See cb() for the two-stage rationale; method='analytic' is the delta-method / generated-regressor correction H^{-1} + sum_i v_i (dphi/dt_i)(dphi/dt_i)'.

mean(Z: Any = None) → float

Mean of the fitted life model.

For an accelerated model the mean life at stress Z is obtained by integrating the survival function (the regression model has no closed mean). For a step-stress model it is the reference-stress mean divided by the acceleration factor at a constant stress, and the integral of the survival function under a StepSchedule.

path(x: ArrayLike, unit: Any) → NDArray

Evaluate the fitted degradation path of unit at x.

For a step-stress (acceleration="clock") model x is calendar time: the path is evaluated on the unit’s reference-stress clock, from its recorded stress history (the last stress held beyond its last measurement).

path_param_link_mean(Z: Any) → NDArray

Mean of the path parameters at stress Z, on the link scale.

For a model fitted with links this is \(D(z)\,\gamma\) – the population mean of the link-scale path parameters eta for a unit tested at stress Z. The parameters are in path order, named like the intercepts in path_param_fixed_names ("log(b)" for a log-linked b). The between-unit covariance around it is path_param_link_cov, the same at every stress.

Parameters:: Z (array like or DataFrame) – One stress row, with as many covariates as the model was fitted with (a one-row DataFrame for a model fitted with fit_from_df). A missing (nan) covariate makes the parameters that depend on it nan.
Returns:: The link-scale mean path parameters at Z.
Return type:: ndarray

path_param_median(Z: Any) → NDArray

Median path parameters at stress Z, on their natural scale.

Each link is monotone and each link-scale parameter is normal, so mapping the link-scale mean through the links gives every parameter’s population median exactly: \(h(D(z)\,\gamma)\). For a log-linked rate that is the geometric-mean rate at Z (the rate’s population mean is larger, by the log-normal factor). Evaluate the typical path at a stress with model.path_model.path(t, *model.path_param_median(Z)).

Parameters:: Z (array like or DataFrame) – One stress row, with as many covariates as the model was fitted with (a one-row DataFrame for a model fitted with fit_from_df). A missing (nan) covariate makes the parameters that depend on it nan.
Returns:: The median path parameters at Z, in path order.
Return type:: ndarray

plot(ax: Any = None) → Any

Plot the degradation data, the fitted per-unit paths (extended to each unit’s pseudo failure time), and the failure threshold.

Parameters:: ax (matplotlib axes, optional) – An axes object to draw the plot on. Creates a new one if not provided.
Returns:: An axes object with the plot.
Return type:: matplotlib axes

predict_failure_time(x: ArrayLike, y: ArrayLike, Z: Any = None, Z_future: Any = None) → float

Estimate the failure time of a new unit from its (partial) degradation trajectory.

Fits this model’s path model to the new unit’s measurements and extrapolates the fitted path to this model’s failure threshold, exactly as was done for each unit during fitting.

Parameters:

x (array like) – Times at which the new unit’s measurements were taken.
y (array like) – The new unit’s degradation measurements.
Z (array like, optional) – For a step-stress (acceleration="clock") model, required: the new unit’s stress history, one row per measurement (the stress over the interval ending at it), or one row for a constant stress. The path is fitted on the unit’s clock.
Z_future (array like or StepSchedule, optional) – For a step-stress model, the stress from the last measurement on: one row, or a StepSchedule whose time zero is the last measurement. Defaults to holding the last stress.

Returns:

The time at which the new unit’s fitted path reaches the threshold. This can be smaller than the last observed time if the trajectory has already crossed the threshold. A trajectory already past the threshold at its first measurement returns the non-positive time at which its fitted path crossed (0 if the path is past the threshold throughout, and for a step-stress model, whose clock starts at zero). Returns nan (with a warning) if the fitted path never reaches the threshold.

Return type:

float

predict_remaining_life(x: ArrayLike, y: ArrayLike, Z: Any = None, Z_future: Any = None) → float

Estimate the remaining life of a new unit from its (partial) degradation trajectory.

This is predict_failure_time() minus the new unit’s last observed time. A negative value means the fitted path crossed the threshold before the last observation (the unit is predicted to have already failed); nan (with a warning) means the fitted path never reaches the threshold. Z and Z_future are as for predict_failure_time().

predict_rul(x: ArrayLike, y: ArrayLike, *, Z: Any = None, Z_future: Any = None, alpha_ci: float = 0.05, n_samples: int = 10000, random_state: int | None = None) → RULPrediction

Bayesian remaining-useful-life prediction for a new unit.

The population distribution of path parameters estimated at fit time (path_param_mean, path_param_cov) is used as a prior, the new unit’s measurements as the likelihood (with the pooled measurement_var as the noise variance), and the Gaussian posterior of the unit’s path parameters is pushed through the threshold crossing by Monte Carlo. The posterior is exact (conjugate) for path models that are linear in their parameters, and an iterated-linearisation (Laplace) approximation otherwise.

Compared to predict_failure_time(), this shrinks short or noisy trajectories toward the population instead of trusting the raw extrapolation, works from a single measurement, and returns credible intervals. With many measurements the posterior concentrates on the least-squares fit and the two agree.

Parameters:

x (array like) – Times at which the new unit’s measurements were taken. One or more measurements are required.
y (array like) – The new unit’s degradation measurements.
Z (array like, optional) – The stress the new unit runs at. Required for a model whose path parameters were modelled against stress (fitted with links): the prior is then the stress-conditional population, eta ~ N(D(z) gamma, Sigma) on the link scale, rather than the pooled population that mixes the stress levels. The posterior is taken on the link scale (so a log-linked rate stays positive) and pushed through the threshold crossing in the same way. Refused for a model without links. For a step-stress (acceleration="clock") model, required: the new unit’s stress history, one row per measurement (the stress over the interval ending at it, as at fit), or one row for a constant stress. The posterior is then taken on the unit’s clock against the reference-stress population, and each sampled reference-stress failure time is mapped back to calendar time along the unit’s history and Z_future.
Z_future (array like or StepSchedule, optional) – For a step-stress model, the stress from the last measurement on: one row, or a StepSchedule whose time zero is the last measurement. Defaults to holding the last stress. Refused for other models.
alpha_ci (float, optional) – Significance level for the equal-tailed credible intervals, between 0 and 1. Defaults to 0.05 (95% intervals).
n_samples (int, optional) – Number of Monte Carlo posterior samples (at least one). Defaults to 10,000.
random_state (int or numpy.random.Generator, optional) – Seed or generator for reproducible sampling. None (the default) seeds from numpy’s global RNG, so np.random.seed controls it.

Returns:

Posterior medians, credible intervals, failure probabilities, and the parameter posterior.

Return type:

RULPrediction

Examples

Eight units with their own start and rate, then a new unit seen three times:

>>> import numpy as np
>>> from surpyval.degradation import DegradationAnalysis
>>> rng = np.random.default_rng(1)
>>> x = np.tile(np.arange(100.0, 1100.0, 100.0), 8)
>>> i = np.repeat(np.arange(8), 10)
>>> a = np.repeat(rng.normal(10.0, 3.0, 8), 10)
>>> b = np.repeat(rng.normal(0.3, 0.05, 8), 10)
>>> y = a + b * x + rng.normal(0, 3.0, x.size)
>>> model = DegradationAnalysis.fit(x, y, i, threshold=450)
>>> pred = model.predict_rul(
...     [100.0, 200.0, 300.0], [42.0, 71.0, 99.0], random_state=0
... )
>>> round(pred.failure_time), round(pred.rul)
(1473, 1173)
>>> [round(v) for v in pred.failure_time_interval]
[1395, 1562]

qf(p: ArrayLike, Z: Any = None) → NDArray

Quantile function of the fitted life model.

Plain life models expose their own qf; accelerated regression models do not, so the quantile at stress Z is obtained by numerically inverting the survival function, pairing each p with a row of Z as sf() pairs each x (a single row, or a single p, is broadcast). For a step-stress model it is the calendar time at which the clock of Z reaches the reference-stress quantile, \(\tau^{-1}(F_0^{-1}(p))\). A missing (nan) probability or covariate gives nan.

random(size: int, Z: Any = None, random_state: int | None = None) → NDArray

Random pseudo failure times from the fitted life model.

For an accelerated model, size samples are drawn at stress Z by inverse-transform sampling of the fitted survival function (the regression models do not all expose random directly). For a step-stress model the reference-stress quantiles are carried to calendar time along the clock of Z.

Parameters:

size (int) – Number of draws.
Z (array like or StepSchedule, optional) – The stress, as for sf(); required for an accelerated or step-stress model, refused otherwise.
random_state (int or numpy.random.Generator, optional) – Seed or generator for reproducible draws. None (the default) seeds from numpy’s global RNG, so np.random.seed controls it.

sf(x: ArrayLike, Z: Any = None) → NDArray

Survival function of the fitted life model.

For an accelerated-degradation model (fitted with covariates) the stress vector Z at which to evaluate life is required. For a step-stress model (acceleration="clock") Z is one stress row or a StepSchedule stress profile, and life is the reference-stress life at the clock time, S(t) = S0(tau(t)); the same holds for every life method below.

to_dict() → dict

Serialise this fitted degradation model to a plain, JSON-serialisable dict.

Everything needed to rebuild the model is stored: the raw measurement data, the path model (by name) and its per-unit fitted parameters, the pseudo failure times and censor flags, the fitted life model (its own to_dict), and the population summaries. The restored model reproduces the life predictions and per-unit paths, and (because the raw data is kept) its bootstrap confidence bounds too.

See also

from_dict, to_json, from_json

to_json(fp: str | PathLike | None = None, with_data: bool = False) → str | None

Write to_dict() to fp as strict JSON, or return it.

Parameters:

fp (str or os.PathLike, optional) – The file to write. Without it the JSON is returned as a string (as pandas.DataFrame.to_json does), which from_json also reads.
with_data (bool, optional) – Write to_dict(with_data=True), which also stores the fitted data, for the models whose to_dict takes with_data (the univariate Parametric and NonParametric); a TypeError for any other model. Defaults to False.

class surpyval.degradation.degradation_analysis.RULPrediction(failure_time: float, failure_time_interval: tuple[float, float], rul: float, rul_interval: tuple[float, float], prob_failed: float, prob_never_fails: float, posterior_mean: NDArray, posterior_cov: NDArray, alpha_ci: float, samples: NDArray)

Bases: object

Posterior failure-time / remaining-useful-life prediction for a new unit, returned by DegradationModel.predict_rul().

All summaries come from Monte Carlo samples of the new unit’s path parameters drawn from their Gaussian posterior and pushed through the path model’s threshold crossing. Samples whose path never reaches the threshold contribute inf failure times, so the median and interval endpoints are inf when that much of the posterior mass never fails. Samples whose path is already past the threshold at the unit’s first measurement (it crossed at or before time zero) have failed: they contribute a failure time of 0, so a trajectory that starts past the threshold has prob_failed = 1, failure_time = 0 and a remaining life of minus its age.

Parameters:

failure_time (float) – Posterior median of the unit’s failure time (measured from the unit’s time zero, like the fitted life model).
failure_time_interval (tuple of float) – Equal-tailed 1 - alpha_ci credible interval for the failure time.
rul (float) – Posterior median remaining useful life: failure time minus the unit’s last observed time. Negative means the unit has most likely already crossed the threshold.
rul_interval (tuple of float) – Equal-tailed 1 - alpha_ci credible interval for the remaining useful life.
prob_failed (float) – Posterior probability that the unit’s path has already crossed the threshold (failure time at or before its last observed time).
prob_never_fails (float) – Posterior probability that the unit’s path never reaches the threshold (a path already past it has failed, not “never fails”).
posterior_mean (ndarray) – The Gaussian posterior of the unit’s path parameters. For a model whose path parameters were modelled against stress (links) these are on the link scale, in the order of the model’s path_param_fixed_names intercepts ("log(b)" for a log-linked b); otherwise on the natural scale.
posterior_cov (ndarray) – The Gaussian posterior of the unit’s path parameters. For a model whose path parameters were modelled against stress (links) these are on the link scale, in the order of the model’s path_param_fixed_names intercepts ("log(b)" for a log-linked b); otherwise on the natural scale.
alpha_ci (float) – The interval significance level used.
samples (ndarray) – The Monte Carlo failure-time samples (inf where the sampled path never reaches the threshold, 0 where it is already past the threshold at the first measurement).

Examples

A new unit, measured three times, of a population of eight fitted units:

>>> import numpy as np
>>> from surpyval.degradation import DegradationAnalysis
>>> rng = np.random.default_rng(1)
>>> x = np.tile(np.arange(100.0, 1100.0, 100.0), 8)
>>> i = np.repeat(np.arange(8), 10)
>>> a = np.repeat(rng.normal(10.0, 3.0, 8), 10)
>>> b = np.repeat(rng.normal(0.3, 0.05, 8), 10)
>>> y = a + b * x + rng.normal(0, 3.0, x.size)
>>> model = DegradationAnalysis.fit(x, y, i, threshold=450)
>>> pred = model.predict_rul(
...     [100.0, 200.0, 300.0], [42.0, 71.0, 99.0], random_state=0
... )
>>> round(pred.rul), [round(v) for v in pred.rul_interval]
(1173, [1095, 1262])
>>> pred.prob_failed
0.0

class surpyval.degradation.degradation_analysis.InducedFailureDistribution(samples: NDArray, threshold: float, path_name: str, stress: list[float] | None = None)

Bases: SerialisableMixin

The population failure-time distribution induced by the degradation path model – the Lu-Meeker approach.

Where the fitted life_model fits a lifetime distribution to each unit’s (noisy) extrapolated pseudo failure time, this instead derives the population life directly from the fitted path-parameter distribution: path parameters are drawn theta ~ N(path_param_mean, path_param_cov) and each draw is pushed through the path model’s inv_path(threshold) to a failure time by Monte Carlo. It is produced by DegradationModel.induced_life().

Draws whose path never reaches the threshold are recorded as inf – a defective (“never fails”) mass exposed as prob_never_fails – so the quantiles and the mean are inf once they reach into that mass. Draws whose path is already past the threshold at the earliest measurement time (it crossed at or before time zero) have failed from the start and are recorded as 0, an atom of failures at time zero.

Use it as a diagnostic: overlay induced.ff(t) on the model’s own ff(t) (the pseudo-failure fit); close agreement is evidence that the path model and its population summary are consistent with the pseudo-failure lifetime fit.

Parameters:

samples (numpy array) – The Monte-Carlo failure-time draws (inf where the path never reaches the threshold, 0 where it is past it from the start).
threshold (float) – The degradation failure threshold used.
path_name (str) – Name of the degradation path model.
stress (list of float, optional) – The stress row the distribution was induced at, for a model whose path parameters depend on stress; None for the plain population.

Examples

The induced life of a fitted population of eight units, next to the Weibull fitted to their pseudo failure times:

>>> import numpy as np
>>> from surpyval.degradation import DegradationAnalysis
>>> rng = np.random.default_rng(1)
>>> x = np.tile(np.arange(100.0, 1100.0, 100.0), 8)
>>> i = np.repeat(np.arange(8), 10)
>>> a = np.repeat(rng.normal(10.0, 3.0, 8), 10)
>>> b = np.repeat(rng.normal(0.3, 0.05, 8), 10)
>>> y = a + b * x + rng.normal(0, 3.0, x.size)
>>> model = DegradationAnalysis.fit(x, y, i, threshold=450)
>>> induced = model.induced_life(random_state=0)
>>> induced
InducedFailureDistribution(Linear path, threshold=450, median=1451.95,
prob_never_fails=0)
>>> induced.sf([1200, 1500]).round(4)
array([0.9949, 0.3489])
>>> model.sf([1200, 1500]).round(4)
array([0.9505, 0.4147])

ff(x: ArrayLike) → NDArray: Failure probability P(T <= x) from the Monte-Carlo draws (nan at a missing time).

classmethod from_dict(model_dict: dict) → InducedFailureDistribution: Rebuild an induced failure-time distribution from a dict.

classmethod from_json(fp: str | PathLike) → Any: Load a model from a JSON file written by to_json(), or from the JSON text it returned (a string starting with {).

mean() → float: Mean failure time (inf if any draw never fails).

median() → float: Median failure time.

qf(p: ArrayLike) → NDArray: Quantile of the induced distribution (inf in the never-fails mass, nan for a missing probability).

random(size: int, random_state: int | None = None) → NDArray: Draw failure times by resampling the Monte-Carlo population.

sf(x: ArrayLike) → NDArray: Survival function P(T > x).

to_dict() → dict

Serialise this induced failure-time distribution to a plain dict.

Stores the Monte-Carlo samples (the inf never-fails draws are written as null so the result is valid JSON), the threshold and the path model’s name.

to_json(fp: str | PathLike | None = None, with_data: bool = False) → str | None

Write to_dict() to fp as strict JSON, or return it.

Parameters:

fp (str or os.PathLike, optional) – The file to write. Without it the JSON is returned as a string (as pandas.DataFrame.to_json does), which from_json also reads.
with_data (bool, optional) – Write to_dict(with_data=True), which also stores the fitted data, for the models whose to_dict takes with_data (the univariate Parametric and NonParametric); a TypeError for any other model. Defaults to False.

Path Models

The shape of each unit’s degradation over time. Pass one to DegradationAnalysis.fit by name (path="linear") or as an instance (path=LinearPath); the instances LinearPath, QuadraticPath, ExponentialPath, … are of the classes below, and PATH_MODELS maps each accepted name ("linear", "quadratic", "exponential", "offset-exponential", "power", "logarithmic", "lloyd-lipow", "gompertz", "michaelis-menten") to its instance. Subclass PathModel for a new shape.

class surpyval.degradation.path_models.PathModel

Bases: ABC

Base class for degradation path models.

A path model is a deterministic function of time with a small number of parameters that is fitted, per unit, to that unit’s degradation measurements. Subclass this (implementing path, inv_path and the name/parameter_names attributes, and either a _initial_guess(x, y) starting point for the default least-squares fit or fit itself) to use a custom degradation path with DegradationAnalysis. A subclass that still names its parameters param_names (before v0.22) works until v0.23, with a DeprecationWarning.

check_data(x: NDArray, y: NDArray) → None: Raise ValueError if the data is outside the model domain.

fit(x: ArrayLike, y: ArrayLike) → NDArray

Fit the path parameters to one unit’s measurements by (nonlinear) least squares.

Parameters:

x (array_like) – The unit’s measurement times.
y (array_like) – Its degradation measurements.

Returns:

The fitted parameters, in the order of parameter_names.

Return type:

numpy array

Examples

>>> from surpyval.degradation import LinearPath
>>> LinearPath.fit([1, 2, 3, 4], [10.5, 12.1, 13.4, 15.2]).round(4)
array([8.95, 1.54])

abstractmethod inv_path(y: ArrayLike, *params: float) → NDArray

Time at which the path reaches level y.

Returns a non-finite value (nan or inf) or a non-positive value when the path never reaches y at a positive finite time.

jacobian(x: ArrayLike, *params: float) → NDArray

Partial derivatives of path with respect to the parameters, evaluated at time(s) x: a (len(x), n_params) matrix.

Used to estimate the least-squares estimation covariance of per-unit fitted parameters. The base implementation uses central finite differences; the built-in models override it with the analytic derivatives.

linear_in_parameters: bool = False: True when path is linear in its parameters, i.e. path(x, *theta) == jacobian(x) @ theta with a Jacobian that does not depend on theta. Enables exact conjugate posterior updates and REML population estimation.

abstractmethod path(x: ArrayLike, *params: float) → NDArray: Evaluate the degradation path at time(s) x.

class surpyval.degradation.path_models.LinearPath_

Bases: PathModel

Linear degradation path: y = a + b * x.

check_data(x: NDArray, y: NDArray) → None: Raise ValueError if the data is outside the model domain.

fit(x: ArrayLike, y: ArrayLike) → NDArray

Fit the path parameters to one unit’s measurements by (nonlinear) least squares.

Parameters:

x (array_like) – The unit’s measurement times.
y (array_like) – Its degradation measurements.

Returns:

The fitted parameters, in the order of parameter_names.

Return type:

numpy array

Examples

>>> from surpyval.degradation import LinearPath
>>> LinearPath.fit([1, 2, 3, 4], [10.5, 12.1, 13.4, 15.2]).round(4)
array([8.95, 1.54])

inv_path(y: ArrayLike, *params: float) → NDArray

Time at which the path reaches level y.

Returns a non-finite value (nan or inf) or a non-positive value when the path never reaches y at a positive finite time.

jacobian(x: ArrayLike, *params: float) → NDArray

Partial derivatives of path with respect to the parameters, evaluated at time(s) x: a (len(x), n_params) matrix.

Used to estimate the least-squares estimation covariance of per-unit fitted parameters. The base implementation uses central finite differences; the built-in models override it with the analytic derivatives.

linear_in_parameters: bool = True: True when path is linear in its parameters, i.e. path(x, *theta) == jacobian(x) @ theta with a Jacobian that does not depend on theta. Enables exact conjugate posterior updates and REML population estimation.

path(x: ArrayLike, *params: float) → NDArray: Evaluate the degradation path at time(s) x.

class surpyval.degradation.path_models.QuadraticPath_

Bases: PathModel

Quadratic degradation path: y = a + b * x + c * x**2.

check_data(x: NDArray, y: NDArray) → None: Raise ValueError if the data is outside the model domain.

fit(x: ArrayLike, y: ArrayLike) → NDArray

Fit the path parameters to one unit’s measurements by (nonlinear) least squares.

Parameters:

x (array_like) – The unit’s measurement times.
y (array_like) – Its degradation measurements.

Returns:

The fitted parameters, in the order of parameter_names.

Return type:

numpy array

Examples

>>> from surpyval.degradation import LinearPath
>>> LinearPath.fit([1, 2, 3, 4], [10.5, 12.1, 13.4, 15.2]).round(4)
array([8.95, 1.54])

inv_path(y: ArrayLike, *params: float) → NDArray: First positive time at which the parabola reaches y.

jacobian(x: ArrayLike, *params: float) → NDArray

Partial derivatives of path with respect to the parameters, evaluated at time(s) x: a (len(x), n_params) matrix.

Used to estimate the least-squares estimation covariance of per-unit fitted parameters. The base implementation uses central finite differences; the built-in models override it with the analytic derivatives.

linear_in_parameters: bool = True: True when path is linear in its parameters, i.e. path(x, *theta) == jacobian(x) @ theta with a Jacobian that does not depend on theta. Enables exact conjugate posterior updates and REML population estimation.

path(x: ArrayLike, *params: float) → NDArray: Evaluate the degradation path at time(s) x.

class surpyval.degradation.path_models.ExponentialPath_

Bases: PathModel

Exponential degradation path: y = a * exp(b * x).

check_data(x: NDArray, y: NDArray) → None: Raise ValueError if the data is outside the model domain.

fit(x: ArrayLike, y: ArrayLike) → NDArray

Fit the path parameters to one unit’s measurements by (nonlinear) least squares.

Parameters:

x (array_like) – The unit’s measurement times.
y (array_like) – Its degradation measurements.

Returns:

The fitted parameters, in the order of parameter_names.

Return type:

numpy array

Examples

>>> from surpyval.degradation import LinearPath
>>> LinearPath.fit([1, 2, 3, 4], [10.5, 12.1, 13.4, 15.2]).round(4)
array([8.95, 1.54])

inv_path(y: ArrayLike, *params: float) → NDArray

Time at which the path reaches level y.

Returns a non-finite value (nan or inf) or a non-positive value when the path never reaches y at a positive finite time.

jacobian(x: ArrayLike, *params: float) → NDArray

Partial derivatives of path with respect to the parameters, evaluated at time(s) x: a (len(x), n_params) matrix.

Used to estimate the least-squares estimation covariance of per-unit fitted parameters. The base implementation uses central finite differences; the built-in models override it with the analytic derivatives.

linear_in_parameters: bool = False: True when path is linear in its parameters, i.e. path(x, *theta) == jacobian(x) @ theta with a Jacobian that does not depend on theta. Enables exact conjugate posterior updates and REML population estimation.

path(x: ArrayLike, *params: float) → NDArray: Evaluate the degradation path at time(s) x.

class surpyval.degradation.path_models.OffsetExponentialPath_

Bases: PathModel

Offset exponential degradation path: y = a + b * exp(c * x).

Covers exponential growth or decay toward/away from the asymptote a; with a = 0 it reduces to the exponential path.

check_data(x: NDArray, y: NDArray) → None: Raise ValueError if the data is outside the model domain.

fit(x: ArrayLike, y: ArrayLike) → NDArray

Fit the path parameters to one unit’s measurements by (nonlinear) least squares.

Parameters:

x (array_like) – The unit’s measurement times.
y (array_like) – Its degradation measurements.

Returns:

The fitted parameters, in the order of parameter_names.

Return type:

numpy array

Examples

>>> from surpyval.degradation import LinearPath
>>> LinearPath.fit([1, 2, 3, 4], [10.5, 12.1, 13.4, 15.2]).round(4)
array([8.95, 1.54])

inv_path(y: ArrayLike, *params: float) → NDArray

Time at which the path reaches level y.

Returns a non-finite value (nan or inf) or a non-positive value when the path never reaches y at a positive finite time.

jacobian(x: ArrayLike, *params: float) → NDArray

Partial derivatives of path with respect to the parameters, evaluated at time(s) x: a (len(x), n_params) matrix.

Used to estimate the least-squares estimation covariance of per-unit fitted parameters. The base implementation uses central finite differences; the built-in models override it with the analytic derivatives.

linear_in_parameters: bool = False: True when path is linear in its parameters, i.e. path(x, *theta) == jacobian(x) @ theta with a Jacobian that does not depend on theta. Enables exact conjugate posterior updates and REML population estimation.

path(x: ArrayLike, *params: float) → NDArray: Evaluate the degradation path at time(s) x.

class surpyval.degradation.path_models.PowerPath_

Bases: PathModel

Power degradation path: y = a * x**b.

check_data(x: NDArray, y: NDArray) → None: Raise ValueError if the data is outside the model domain.

fit(x: ArrayLike, y: ArrayLike) → NDArray

Fit the path parameters to one unit’s measurements by (nonlinear) least squares.

Parameters:

x (array_like) – The unit’s measurement times.
y (array_like) – Its degradation measurements.

Returns:

The fitted parameters, in the order of parameter_names.

Return type:

numpy array

Examples

>>> from surpyval.degradation import LinearPath
>>> LinearPath.fit([1, 2, 3, 4], [10.5, 12.1, 13.4, 15.2]).round(4)
array([8.95, 1.54])

inv_path(y: ArrayLike, *params: float) → NDArray

Time at which the path reaches level y.

Returns a non-finite value (nan or inf) or a non-positive value when the path never reaches y at a positive finite time.

jacobian(x: ArrayLike, *params: float) → NDArray

Partial derivatives of path with respect to the parameters, evaluated at time(s) x: a (len(x), n_params) matrix.

Used to estimate the least-squares estimation covariance of per-unit fitted parameters. The base implementation uses central finite differences; the built-in models override it with the analytic derivatives.

linear_in_parameters: bool = False: True when path is linear in its parameters, i.e. path(x, *theta) == jacobian(x) @ theta with a Jacobian that does not depend on theta. Enables exact conjugate posterior updates and REML population estimation.

path(x: ArrayLike, *params: float) → NDArray: Evaluate the degradation path at time(s) x.

class surpyval.degradation.path_models.LogarithmicPath_

Bases: PathModel

Logarithmic degradation path: y = a + b * ln(x).

check_data(x: NDArray, y: NDArray) → None: Raise ValueError if the data is outside the model domain.

fit(x: ArrayLike, y: ArrayLike) → NDArray

Fit the path parameters to one unit’s measurements by (nonlinear) least squares.

Parameters:

x (array_like) – The unit’s measurement times.
y (array_like) – Its degradation measurements.

Returns:

The fitted parameters, in the order of parameter_names.

Return type:

numpy array

Examples

>>> from surpyval.degradation import LinearPath
>>> LinearPath.fit([1, 2, 3, 4], [10.5, 12.1, 13.4, 15.2]).round(4)
array([8.95, 1.54])

inv_path(y: ArrayLike, *params: float) → NDArray

Time at which the path reaches level y.

Returns a non-finite value (nan or inf) or a non-positive value when the path never reaches y at a positive finite time.

jacobian(x: ArrayLike, *params: float) → NDArray

Partial derivatives of path with respect to the parameters, evaluated at time(s) x: a (len(x), n_params) matrix.

Used to estimate the least-squares estimation covariance of per-unit fitted parameters. The base implementation uses central finite differences; the built-in models override it with the analytic derivatives.

linear_in_parameters: bool = True: True when path is linear in its parameters, i.e. path(x, *theta) == jacobian(x) @ theta with a Jacobian that does not depend on theta. Enables exact conjugate posterior updates and REML population estimation.

path(x: ArrayLike, *params: float) → NDArray: Evaluate the degradation path at time(s) x.

class surpyval.degradation.path_models.LloydLipowPath_

Bases: PathModel

Lloyd-Lipow degradation path: y = a - b / x.

check_data(x: NDArray, y: NDArray) → None: Raise ValueError if the data is outside the model domain.

fit(x: ArrayLike, y: ArrayLike) → NDArray

Fit the path parameters to one unit’s measurements by (nonlinear) least squares.

Parameters:

x (array_like) – The unit’s measurement times.
y (array_like) – Its degradation measurements.

Returns:

The fitted parameters, in the order of parameter_names.

Return type:

numpy array

Examples

>>> from surpyval.degradation import LinearPath
>>> LinearPath.fit([1, 2, 3, 4], [10.5, 12.1, 13.4, 15.2]).round(4)
array([8.95, 1.54])

inv_path(y: ArrayLike, *params: float) → NDArray

Time at which the path reaches level y.

Returns a non-finite value (nan or inf) or a non-positive value when the path never reaches y at a positive finite time.

jacobian(x: ArrayLike, *params: float) → NDArray

Partial derivatives of path with respect to the parameters, evaluated at time(s) x: a (len(x), n_params) matrix.

Used to estimate the least-squares estimation covariance of per-unit fitted parameters. The base implementation uses central finite differences; the built-in models override it with the analytic derivatives.

linear_in_parameters: bool = True: True when path is linear in its parameters, i.e. path(x, *theta) == jacobian(x) @ theta with a Jacobian that does not depend on theta. Enables exact conjugate posterior updates and REML population estimation.

path(x: ArrayLike, *params: float) → NDArray: Evaluate the degradation path at time(s) x.

class surpyval.degradation.path_models.GompertzPath_

Bases: PathModel

Gompertz degradation path: y = a * exp(-b * exp(-c * x)).

An S-shaped path approaching the asymptote a.

check_data(x: NDArray, y: NDArray) → None: Raise ValueError if the data is outside the model domain.

fit(x: ArrayLike, y: ArrayLike) → NDArray

Fit the path parameters to one unit’s measurements by (nonlinear) least squares.

Parameters:

x (array_like) – The unit’s measurement times.
y (array_like) – Its degradation measurements.

Returns:

The fitted parameters, in the order of parameter_names.

Return type:

numpy array

Examples

>>> from surpyval.degradation import LinearPath
>>> LinearPath.fit([1, 2, 3, 4], [10.5, 12.1, 13.4, 15.2]).round(4)
array([8.95, 1.54])

inv_path(y: ArrayLike, *params: float) → NDArray

Time at which the path reaches level y.

Returns a non-finite value (nan or inf) or a non-positive value when the path never reaches y at a positive finite time.

jacobian(x: ArrayLike, *params: float) → NDArray

Partial derivatives of path with respect to the parameters, evaluated at time(s) x: a (len(x), n_params) matrix.

Used to estimate the least-squares estimation covariance of per-unit fitted parameters. The base implementation uses central finite differences; the built-in models override it with the analytic derivatives.

linear_in_parameters: bool = False: True when path is linear in its parameters, i.e. path(x, *theta) == jacobian(x) @ theta with a Jacobian that does not depend on theta. Enables exact conjugate posterior updates and REML population estimation.

path(x: ArrayLike, *params: float) → NDArray: Evaluate the degradation path at time(s) x.

class surpyval.degradation.path_models.MichaelisMentenPath_

Bases: PathModel

Michaelis-Menten degradation path: y = a * x / (b + x).

A saturating path rising from zero toward the asymptote a, reaching half of it at x = b.

check_data(x: NDArray, y: NDArray) → None: Raise ValueError if the data is outside the model domain.

fit(x: ArrayLike, y: ArrayLike) → NDArray

Fit the path parameters to one unit’s measurements by (nonlinear) least squares.

Parameters:

x (array_like) – The unit’s measurement times.
y (array_like) – Its degradation measurements.

Returns:

The fitted parameters, in the order of parameter_names.

Return type:

numpy array

Examples

>>> from surpyval.degradation import LinearPath
>>> LinearPath.fit([1, 2, 3, 4], [10.5, 12.1, 13.4, 15.2]).round(4)
array([8.95, 1.54])

inv_path(y: ArrayLike, *params: float) → NDArray

Time at which the path reaches level y.

Returns a non-finite value (nan or inf) or a non-positive value when the path never reaches y at a positive finite time.

jacobian(x: ArrayLike, *params: float) → NDArray

Partial derivatives of path with respect to the parameters, evaluated at time(s) x: a (len(x), n_params) matrix.

Used to estimate the least-squares estimation covariance of per-unit fitted parameters. The base implementation uses central finite differences; the built-in models override it with the analytic derivatives.

linear_in_parameters: bool = False: True when path is linear in its parameters, i.e. path(x, *theta) == jacobian(x) @ theta with a Jacobian that does not depend on theta. Enables exact conjugate posterior updates and REML population estimation.

path(x: ArrayLike, *params: float) → NDArray: Evaluate the degradation path at time(s) x.

surpyval.degradation.path_models.get_path_model(path: str | PathModel) → PathModel

Resolve path to a PathModel instance.

Accepts a PathModel instance (returned unchanged) or one of the registered names in PATH_MODELS (case-insensitive): "linear", "quadratic", "exponential", "offset-exponential", "power", "logarithmic", "lloyd-lipow", "gompertz", "michaelis-menten". A built-in model’s display name (e.g. "Offset Exponential") is accepted too. ("best" — automatic selection — is handled by DegradationAnalysis.fit, not here.)

Examples

>>> import numpy as np
>>> from surpyval.degradation import get_path_model
>>> power = get_path_model("power")
>>> power.name, power.parameter_names
('Power', ['a', 'b'])
>>> power.path(np.array([1.0, 4.0]), 2.0, 0.5)
array([2., 4.])
>>> get_path_model("Offset Exponential").name
'Offset Exponential'

Stress-Dependent Path Parameters

For accelerated degradation tests whose mechanism depends on stress (links in DegradationAnalysis.fit): the path parameters are modelled on a link scale, eta_i = D(z_i) gamma + u_i, so a log-linked rate with Z = 1/T follows the Arrhenius relationship.

class surpyval.degradation.stress.LinkedPathModel(base: PathModel, links: dict[str, str])

Bases: PathModel

A path model reparameterised onto a link scale.

Wraps base so that its parameters theta are replaced by eta with theta = h(eta) elementwise: h is the identity for most parameters and exp for those given a "log" link. path, inv_path, jacobian and fit all take and return link-scale parameters, so the wrapped model is a drop-in path model for the per-unit fits and the population (REML) machinery, which then estimate the population of eta.

Parameters:

base (PathModel) – The path model on its natural scale.
links (dict) – {parameter name: "identity" | "log"} for the parameters whose link is not the identity (a parameter left out gets the identity link).

Examples

The linear path a + b t with its slope on a log link, so that b stays positive:

>>> import numpy as np
>>> from surpyval.degradation import LinkedPathModel, get_path_model
>>> linked = LinkedPathModel(get_path_model("linear"), {"b": "log"})
>>> linked.parameter_names
['a', 'log(b)']
>>> eta = linked.to_link([10.0, 0.3])
>>> eta.round(4)
array([10.   , -1.204])
>>> linked.path(np.array([0.0, 10.0]), *eta)
array([10., 13.])

check_data(x: NDArray, y: NDArray) → None: Raise ValueError if the data is outside the model domain.

fit(x: ArrayLike, y: ArrayLike) → NDArray: The base model’s per-unit fit, returned on the link scale.

inv_path(y: ArrayLike, *params: float) → NDArray

Time at which the path reaches level y.

Returns a non-finite value (nan or inf) or a non-positive value when the path never reaches y at a positive finite time.

jacobian(x: ArrayLike, *params: float) → NDArray

Partial derivatives of path with respect to the parameters, evaluated at time(s) x: a (len(x), n_params) matrix.

Used to estimate the least-squares estimation covariance of per-unit fitted parameters. The base implementation uses central finite differences; the built-in models override it with the analytic derivatives.

linear_in_parameters: bool = False: True when path is linear in its parameters, i.e. path(x, *theta) == jacobian(x) @ theta with a Jacobian that does not depend on theta. Enables exact conjugate posterior updates and REML population estimation.

path(x: ArrayLike, *params: float) → NDArray: Evaluate the degradation path at time(s) x.

to_link(theta: ArrayLike) → NDArray: Link-scale parameters eta = h^-1(theta); the parameters run along the last axis, as for to_natural().

to_natural(eta: ArrayLike) → NDArray

Natural-scale parameters theta = h(eta).

eta is one parameter vector, or an array of them with the parameters along the last axis (one row per Monte-Carlo draw).

surpyval.degradation.stress.stress_design(z: ArrayLike, links: dict[str, str], parameter_names: list[str]) → NDArray

The fixed-effects design D(z) of one unit: a (p, m) matrix mapping the fixed effects gamma to the unit’s link-scale path parameter means, eta_mean = D(z) gamma.

The columns are laid out parameter by parameter, in path order: an intercept column for every parameter, followed – for the stress-dependent parameters named in links – by one column per covariate carrying that unit’s stress values. So m = p + q * len(links) with q covariates.

surpyval.degradation.stress.fixed_effect_names(linked_parameter_names: list[str], parameter_names: list[str], links: dict[str, str], n_cov: int) → list[str]: Labels for gamma matching stress_design()’s columns: the link-scale parameter name for each intercept and "<name>:Z<j>" for each stress coefficient.

Step-Stress: the Accelerated Clock

For tests whose stress changes during a unit’s test (acceleration="clock" in DegradationAnalysis.fit): stress speeds up the clock of every unit’s path, AF(z) = exp(gamma' (z - stress_ref)), and the path is the ordinary path model on the reference-stress time the unit has aged. The fitted DegradationModel then carries gamma and stress_ref, its life methods take the stress as one row or a StepSchedule, and its trajectory methods take the unit’s stress history Z and a planned Z_future. How the stress coefficients are estimated:

Step-stress general-path degradation: estimating the accelerated clock.

With acceleration="clock" in DegradationAnalysis.fit() stress speeds up the clock of every unit’s degradation path, so unit i’s measurements are

y_ij = g(tau_i(t_ij); theta_i) + eps_ij, theta_i ~ MVN(mu, Sigma),

with tau_i(t) = sum_k AF(z_ik) dt_ik the reference-stress time the unit has aged by t and AF(z) = exp(gamma' (z - z_ref)). Z has one row per measurement: the stress applied over the interval that ends at that measurement (the first interval starts at time zero).

Once gamma is known the model is the ordinary general-path model on tau, so everything downstream reuses the existing pipeline. This module estimates gamma:

profile_least_squares() – the two-stage route. For a trial gamma every unit’s path is refitted on its tau and the residual sums of squares are pooled; gamma minimises the total. A unit held at one stress can absorb any acceleration into its own path parameters (every built-in path family is closed under rescaling time), so only units whose stress changes during the test carry information here.
mixed_model_estimate() – the mixed-model route. The units share one population of path parameters, so units at different constant stresses identify gamma too. A Lindstrom-Bates FOCE iteration with gamma among the fixed effects (the path linearised in theta_i and gamma) finds the neighbourhood, and the FOCE-approximate profile likelihood of gamma is then maximised directly.

class surpyval.degradation.step_stress.ClockUnit(x: NDArray, y: NDArray, dt: NDArray, s: NDArray)

Bases: object

One unit’s time-ordered data with the stress over each interval.

tau(g: NDArray) → NDArray: Reference-stress time at each measurement.

tau_gradient(g: NDArray) → NDArray: d tau / d g at each measurement, (n, q).

surpyval.degradation.step_stress.clock_units(x: NDArray, y: NDArray, i: NDArray, units: NDArray, Z: NDArray, z_ref: NDArray, scale: NDArray) → list[ClockUnit]: Split the data by unit, time-ordered, with scaled interval stresses.

surpyval.degradation.step_stress.mixed_model_estimate(units: list[ClockUnit], path_model: Any, g_init: NDArray, theta_init: NDArray, mean_init: NDArray, cov_init: NDArray, sigma2_init: float) → tuple[NDArray, bool, tuple]

Mixed-model estimate of the scaled stress coefficients.

The joint FOCE iteration (g among the fixed effects) moves quickly to the neighbourhood of the estimate, but crawls along the ridge a unit held at one stress creates – its rate and g nearly trade off, and the linearisation follows that curved ridge in small steps. So it is followed by a direct maximisation of the profile likelihood of g: for each trial g the population is refitted by FOCE with g held fixed, and the approximate marginal likelihood compared. g changes the fixed-effects design, so the profile uses the plain likelihood rather than REML, whose adjustment term is not comparable across designs.

Returns (g, converged, population) with population the REML estimate (mu, Sigma, sigma2, converged) of the reference-stress population of path parameters at that g.

surpyval.degradation.step_stress.path_time_derivative(path_model: Any, tau: NDArray, theta: NDArray) → NDArray: d g(tau; theta) / d tau by finite differences (one-sided at 0).

surpyval.degradation.step_stress.profile_least_squares(units: list[ClockUnit], path_model: Any, q: int) → NDArray: Two-stage estimate of the scaled stress coefficients: minimise the pooled per-unit residual sum of squares of the paths refitted on tau(g). For one covariate a grid locates the basin before a bounded scalar search; otherwise a local search starts from g = 0.

Stochastic Process Models

Where a path model treats degradation as a deterministic curve with noise, these treat it as a stochastic process in its own right: the Wiener process for degradation that can go down as well as up, and the gamma process for monotone accumulation such as wear or crack growth. Both give a first-passage distribution to the threshold in closed form, and so a remaining-useful-life prediction with bounds. Both also take a stress Z (one row per measurement), which may change between or during units’ tests: stress accelerates the process clock, and the life under any stress profile stays in closed form.

class surpyval.degradation.process_models.WienerProcess

Bases: object

Fitter for the Wiener-process degradation model (see WienerProcessModel).

classmethod fit(x: ArrayLike, y: ArrayLike, i: ArrayLike, threshold: float, Z: ArrayLike | None = None, stress_ref: ArrayLike | None = None) → WienerProcessModel

Fit a Wiener-process degradation model by maximum likelihood.

Parameters:

x (array_like) – Measurement times.
y (array_like) – Degradation measurements.
i (array_like) – Unit identifier for each measurement.
threshold (float) – The degradation level defining failure.
Z (array_like, optional) – Stress covariates, one row per measurement, for an accelerated or step-stress test. The stress may differ between units and change during a unit’s test: the stress on a measurement is the stress applied since the previous one. Stress speeds up the process clock by exp(gamma' (z - stress_ref)), scaling the drift and the variance per unit time together (the time-scale transformation of Whitmore and Schenkelberg, 1997); gamma is estimated with mu and sigma by maximum likelihood. At least two distinct stress levels are needed.
stress_ref (array_like, optional) – The reference stress at which the fitted mu and sigma apply. Defaults to the mean stress over the measurement intervals; pass the use conditions to read the model at them.

Returns:

The fitted model, whose life-distribution methods (sf, ff, mean, …) give the first-passage time to threshold.

Return type:

WienerProcessModel

Examples

Five units whose degradation drifts upwards at 0.5 per unit time with Brownian noise:

>>> import numpy as np
>>> from surpyval.degradation import WienerProcess
>>> rng = np.random.default_rng(1)
>>> t = np.tile(np.arange(0, 110, 10.0), 5)  # 5 units, 11 readings
>>> i = np.repeat(np.arange(5), 11)
>>> steps = rng.normal(0.5 * 10, 1.0 * np.sqrt(10), size=(5, 10))
>>> y = np.hstack([np.r_[0.0, np.cumsum(s)] for s in steps])
>>> model = WienerProcess.fit(t, y, i, threshold=100)
>>> model
Wiener Process Degradation Model
================================
Drift (mu)          : 0.488591
Diffusion (sigma)   : 0.88113
Threshold           : 100
Mean time to failure: 204.67
>>> model.sf([150, 200]).round(4)
array([0.9922, 0.548 ])

classmethod fit_from_df(df: DataFrame, x_col: str = 'x', y_col: str = 'y', i_col: str = 'i', Z_cols: str | list[str] | None = None, **fit_kwargs: Any) → WienerProcessModel

Fit a Wiener-process degradation model from the columns of a DataFrame.

Parameters:

df (DataFrame) – The degradation data, one row per measurement.
x_col (str, optional) – The columns of the measurement times, the measurements and the unit identifiers. Default "x", "y" and "i". Their v0.21 names x, y and i still work, with a DeprecationWarning, until v0.23 (every DataFrame entry point names its columns with a _col suffix, principle 21).
y_col (str, optional) – The columns of the measurement times, the measurements and the unit identifiers. Default "x", "y" and "i". Their v0.21 names x, y and i still work, with a DeprecationWarning, until v0.23 (every DataFrame entry point names its columns with a _col suffix, principle 21).
i_col (str, optional) – The columns of the measurement times, the measurements and the unit identifiers. Default "x", "y" and "i". Their v0.21 names x, y and i still work, with a DeprecationWarning, until v0.23 (every DataFrame entry point names its columns with a _col suffix, principle 21).
Z_cols (str or list of str, optional) – The stress column(s), passed to fit() as Z. Their names are recorded on the model (as Z_cols, kept by to_dict), so its methods also take the stress as a one-row DataFrame and select these columns by name.
**fit_kwargs – The remaining arguments of fit(): threshold (required), and optionally stress_ref and the others.

Returns:

The fitted model.

Return type:

WienerProcessModel

class surpyval.degradation.process_models.WienerProcessModel(mu: float, sigma: float, threshold: float, gamma: Any = None, stress_ref: Any = None)

Bases: FirstPassageProcessModel

A fitted Wiener-process degradation model, W(t) = mu*t + sigma*B(t).

The first passage of the process to the failure threshold (from an assumed degradation of zero at t = 0) is Inverse-Gaussian distributed with mean threshold / mu and shape threshold**2 / sigma**2, and the failure-time methods below evaluate that distribution. A positive drift mu is required for a proper (non-defective) life distribution.

Parameters:

mu (float) – Fitted drift (mean degradation rate).
sigma (float) – Fitted diffusion (volatility) coefficient.
threshold (float) – The degradation level defining failure.
gamma (array like, optional) – For a model fitted with stress Z: the stress coefficients and the reference stress at which mu and sigma apply. At stress z the process clock runs exp(gamma' (z - stress_ref)) times faster, scaling both the drift and the variance per unit time.
stress_ref (array like, optional) – For a model fitted with stress Z: the stress coefficients and the reference stress at which mu and sigma apply. At stress z the process clock runs exp(gamma' (z - stress_ref)) times faster, scaling both the drift and the variance per unit time.

Examples

WienerProcess.fit returns one; it can also be built from known parameters:

>>> from surpyval.degradation import WienerProcessModel
>>> model = WienerProcessModel(mu=0.5, sigma=1.0, threshold=100)
>>> model
Wiener Process Degradation Model
================================
Drift (mu)          : 0.5
Diffusion (sigma)   : 1
Threshold           : 100
Mean time to failure: 200
>>> model.sf([150, 200, 250]).round(4)
array([0.9759, 0.4719, 0.0489])

Hf(x: ArrayLike, Z: Any = None) → NDArray: Cumulative hazard of the first-passage time.

Z_cols: list[str] | None = None: The covariate columns of a model fitted with fit_from_df, by which a DataFrame Z is read; None for a model fitted from arrays.

acceleration_factor(Z: Any) → float

How much faster the process runs at stress Z than at the reference stress: exp(gamma' (z - stress_ref)).

A life at the reference stress divides by this to give the life at Z, and a unit at Z accumulates degradation this many times faster.

Parameters:: Z (array like or DataFrame) – One stress row (nan for a missing value gives nan). A model fitted with fit_from_df also takes a one-row DataFrame, read by column name.

df(x: ArrayLike, Z: Any = None) → NDArray: Density of the first-passage time. Under a stress path it is the reference-stress density at the clock time tau(x) times the clock’s rate, AF of the stress in force at x.

ff(x: ArrayLike, Z: Any = None) → NDArray

Failure (CDF) of the first-passage time to the threshold.

For a model fitted with stress, Z is required: one stress row for a constant stress, or a StepSchedule for a stress profile; a model fitted with fit_from_df also takes a one-row DataFrame, read by column name. The same applies to every method below. A missing (nan) time or stress gives nan.

classmethod from_dict(model_dict: dict) → FirstPassageProcessModel: Rebuild a fitted process model from a to_dict() dict.

classmethod from_json(fp: str | PathLike) → Any: Load a model from a JSON file written by to_json(), or from the JSON text it returned (a string starting with {).

gamma: NDArray | None: Stress coefficients and the reference stress, for a model fitted with Z; both None otherwise.

hf(x: ArrayLike, Z: Any = None) → NDArray: Hazard function of the first-passage time.

property is_accelerated: bool: True for a model fitted with a stress Z.

mean(Z: Any = None) → float: Mean time to failure. At a constant stress it is the reference-stress mean divided by the acceleration factor; under a stress profile it is the integral of the survival function.

predict_rul(current_degradation: float, *, Z: Any = None, alpha_ci: float = 0.05) → ProcessRUL

Remaining useful life given the current degradation level.

Increments are independent for both processes, so the remaining first passage over the residual distance threshold - current_degradation follows the same law as a fresh process; its median and equal-tailed interval are returned.

Parameters:

current_degradation (float) – The unit’s current degradation level (not nan).
Z (array like or StepSchedule, optional) – For a model fitted with stress: the stress the unit will run at from now on – one row for a constant stress, or a StepSchedule whose time zero is now (or, for a model fitted with fit_from_df, a one-row DataFrame). It describes this one unit, so a missing value is refused.
alpha_ci (float, optional) – Tail probability of the returned interval, between 0 and 1. Default 0.05.

Returns:

The median remaining life and its equal-tailed interval.

Return type:

ProcessRUL

qf(p: ArrayLike, Z: Any = None) → NDArray: Quantile (inverse CDF) of the first-passage time (nan for a missing probability or stress).

random(size: int, Z: Any = None, random_state: int | None = None) → NDArray

Draw first-passage (failure) times from the fitted model.

Parameters:

size (int) – Number of draws.
Z (array like or StepSchedule, optional) – The stress, for a model fitted with Z (required then); each reference-stress draw is carried to calendar time along its clock.
random_state (int or numpy.random.Generator, optional) – Seed or generator for reproducible draws. None (the default) seeds from numpy’s global RNG, so np.random.seed controls it.

sf(x: ArrayLike, Z: Any = None) → NDArray: Survival function of the first-passage time.

to_dict() → dict: Serialise this fitted process model to a plain dict.

to_json(fp: str | PathLike | None = None, with_data: bool = False) → str | None

Write to_dict() to fp as strict JSON, or return it.

Parameters:

fp (str or os.PathLike, optional) – The file to write. Without it the JSON is returned as a string (as pandas.DataFrame.to_json does), which from_json also reads.
with_data (bool, optional) – Write to_dict(with_data=True), which also stores the fitted data, for the models whose to_dict takes with_data (the univariate Parametric and NonParametric); a TypeError for any other model. Defaults to False.

class surpyval.degradation.process_models.GammaProcess

Bases: object

Fitter for the Gamma-process degradation model (see GammaProcessModel).

classmethod fit(x: ArrayLike, y: ArrayLike, i: ArrayLike, threshold: float, Z: ArrayLike | None = None, stress_ref: ArrayLike | None = None, resolution: float | None = None, gauge: float | None = None, rounding: str = 'nearest', exact_start: bool = False, gauge_method: str = 'exact') → GammaProcessModel

Fit a Gamma-process degradation model by maximum likelihood.

The degradation must be monotone increasing (all increments non-negative); an increment that decreases raises an error pointing to the Wiener model for non-monotone signals.

A gamma increment is zero with probability zero, so an increment recorded as exactly zero means the change was below the measurement resolution: it enters the likelihood as censored, P(increment <= resolution). Data without zero increments is fitted by the ordinary likelihood.

That treats every non-zero increment as exact, which is fine when the readings are much finer than the increments. When they come from a coarse gauge – a step comparable to the increments – every increment is rounded, and the fit is biased (alpha can come out at double or half its value, see the example below). Pass the gauge step as gauge to fit the quantised likelihood instead: a reading r then means the true level lies somewhere in the gauge bin around r, and the likelihood of a unit is the probability that its whole path passes through its recorded bins. The increments of a path are not independent once rounded (two consecutive increments share the rounding error of the reading between them), so gauge_method="exact" evaluates that path probability by a forward recursion over each unit’s readings, carrying the distribution of the true level across its bin (on a grid of cells, integrated exactly within each cell). gauge_method="independent" is a cheaper approximation that multiplies the exact probabilities of the single recorded increments, ignoring their dependence; it is nearly unbiased too, but less efficient. The true level at a unit’s first reading is taken to be uniform over its bin, unless exact_start says it is known.

A gauge coarser than the scatter a unit’s degradation builds up over the whole test leaves the readings unable to tell a random path from a straight line: alpha and beta then run off to large values together, while their ratio (the mean rate) and the mean life stay well estimated.

Parameters:

x (array_like) – Measurement times.
y (array_like) – Degradation measurements.
i (array_like) – Unit identifier for each measurement.
threshold (float) – The degradation level defining failure.
Z (array_like, optional) – Stress covariates, one row per measurement, for an accelerated or step-stress test. The stress may differ between units and change during a unit’s test: the stress on a measurement is the stress applied since the previous one. Stress speeds up the process clock, so the shape accrues at alpha * exp(gamma' (z - stress_ref)) per unit time with beta unchanged; gamma is estimated with alpha and beta by maximum likelihood. At least two distinct stress levels are needed.
stress_ref (array_like, optional) – The reference stress at which the fitted alpha applies. Defaults to the mean stress over the measurement intervals.
resolution (float, optional) – The measurement resolution: an increment recorded as zero is taken to be somewhere in [0, resolution]. Defaults to the smallest positive increment in the data (for readings rounded to a grid, the grid step). Only used when some increment is zero, and not with gauge, which models the zeros itself.
gauge (float, optional) – The step of the gauge the readings were rounded to. Given, the fit uses the quantised likelihood described above; the differences between readings of a unit must then be whole multiples of it. The default, None, keeps the likelihood above (exact non-zero increments, censored zeros).
rounding ({"nearest", "floor"}, optional) – How the gauge rounds, with gauge: "nearest" (the default) means a reading r stands for a true level in [r - gauge/2, r + gauge/2); "floor" (a gauge that truncates, or a counter that ticks once a whole step is complete) means [r, r + gauge). Only differences between bins enter the likelihood, so the convention changes the fit only together with exact_start.
exact_start (bool, optional) – With gauge: True if each unit’s first reading is its exact true level (for instance new units at zero wear, not read off the gauge), so that later bins are placed relative to it by rounding. The default, False, treats the first reading as a gauge reading like the rest.
gauge_method ({"exact", "independent"}, optional) – With gauge: "exact" (the default) for the path likelihood, "independent" for the cheaper approximation that treats the rounded increments as independent. The exact method costs several times as much: a few hundred increments take a fraction of a second, ten thousand a few seconds.

Returns:

The fitted model, whose life-distribution methods (sf, ff, mean, …) give the first-passage time to threshold.

Return type:

GammaProcessModel

Warns:

UserWarning – “No finite maximum” when every increment is proportional to its time step (noise-free readings), so the likelihood keeps increasing with alpha: the returned alpha and beta are meaningless. (WienerProcess refuses such data.)

Examples

Five units whose wear accumulates in non-negative gamma-distributed increments (mean 0.5 per unit time):

>>> import numpy as np
>>> from surpyval.degradation import GammaProcess
>>> rng = np.random.default_rng(1)
>>> t = np.tile(np.arange(0, 110, 10.0), 5)  # 5 units, 11 readings
>>> i = np.repeat(np.arange(5), 11)
>>> steps = rng.gamma(shape=2.0 * 10, scale=0.25, size=(5, 10))
>>> y = np.hstack([np.r_[0.0, np.cumsum(s)] for s in steps])
>>> model = GammaProcess.fit(t, y, i, threshold=100)
>>> model
Gamma Process Degradation Model
===============================
Shape rate (alpha)  : 2.61045
Rate (beta)         : 5.45561
Threshold           : 100
Mean time to failure: 209.183
>>> model.sf([150, 200]).round(4)
array([1.    , 0.8477])

The same wear read off a gauge with a step of 2. Rounding turns increments of about 5 into 4s and 6s, so fitted as exact increments they look far more variable than they are and alpha halves; the quantised likelihood puts it back near the 2.61 of the unrounded readings:

>>> y_gauge = np.round(y / 2.0) * 2.0
>>> naive = GammaProcess.fit(t, y_gauge, i, threshold=100)
>>> quantised = GammaProcess.fit(t, y_gauge, i, 100, gauge=2.0)
>>> print(round(naive.alpha, 2), round(quantised.alpha, 2))
1.07 2.16

classmethod fit_from_df(df: DataFrame, x_col: str = 'x', y_col: str = 'y', i_col: str = 'i', Z_cols: str | list[str] | None = None, **fit_kwargs: Any) → GammaProcessModel

Fit a Gamma-process degradation model from the columns of a DataFrame.

Parameters:

df (DataFrame) – The degradation data, one row per measurement.
x_col (str, optional) – The columns of the measurement times, the measurements and the unit identifiers. Default "x", "y" and "i". Their v0.21 names x, y and i still work, with a DeprecationWarning, until v0.23 (every DataFrame entry point names its columns with a _col suffix, principle 21).
y_col (str, optional) – The columns of the measurement times, the measurements and the unit identifiers. Default "x", "y" and "i". Their v0.21 names x, y and i still work, with a DeprecationWarning, until v0.23 (every DataFrame entry point names its columns with a _col suffix, principle 21).
i_col (str, optional) – The columns of the measurement times, the measurements and the unit identifiers. Default "x", "y" and "i". Their v0.21 names x, y and i still work, with a DeprecationWarning, until v0.23 (every DataFrame entry point names its columns with a _col suffix, principle 21).
Z_cols (str or list of str, optional) – The stress column(s), passed to fit() as Z. Their names are recorded on the model (as Z_cols, kept by to_dict), so its methods also take the stress as a one-row DataFrame and select these columns by name.
**fit_kwargs – The remaining arguments of fit(): threshold (required), and optionally stress_ref and the others.

Returns:

The fitted model.

Return type:

GammaProcessModel

class surpyval.degradation.process_models.GammaProcessModel(alpha: float, beta: float, threshold: float, gamma: Any = None, stress_ref: Any = None)

Bases: FirstPassageProcessModel

A fitted Gamma-process degradation model with stationary independent increments: over an interval dt the degradation increment is Gamma(shape = alpha * dt, rate = beta). The path is monotone increasing.

Because the path is monotone, the first passage to the failure threshold has failure CDF P(W(t) >= threshold), evaluated with the regularised upper incomplete gamma function.

Parameters:

alpha (float) – Fitted shape rate (shape accrues as alpha * t).
beta (float) – Fitted rate parameter of the increments.
threshold (float) – The degradation level defining failure.
gamma (array like, optional) – For a model fitted with stress Z: the stress coefficients and the reference stress at which alpha applies. At stress z the shape accrues at alpha * exp(gamma' (z - stress_ref)).
stress_ref (array like, optional) – For a model fitted with stress Z: the stress coefficients and the reference stress at which alpha applies. At stress z the shape accrues at alpha * exp(gamma' (z - stress_ref)).

Examples

GammaProcess.fit returns one; it can also be built from known parameters. Wear accruing at a mean alpha / beta = 0.5 per unit time, with failure at 100:

>>> from surpyval.degradation import GammaProcessModel
>>> model = GammaProcessModel(alpha=2.0, beta=4.0, threshold=100)
>>> model.sf([150, 200, 250]).round(4)
array([1.    , 0.5066, 0.    ])
>>> round(model.mean(), 2)
200.25

Hf(x: ArrayLike, Z: Any = None) → NDArray: Cumulative hazard of the first-passage time.

Z_cols: list[str] | None = None: The covariate columns of a model fitted with fit_from_df, by which a DataFrame Z is read; None for a model fitted from arrays.

acceleration_factor(Z: Any) → float

How much faster the process runs at stress Z than at the reference stress: exp(gamma' (z - stress_ref)).

A life at the reference stress divides by this to give the life at Z, and a unit at Z accumulates degradation this many times faster.

Parameters:: Z (array like or DataFrame) – One stress row (nan for a missing value gives nan). A model fitted with fit_from_df also takes a one-row DataFrame, read by column name.

df(x: ArrayLike, Z: Any = None) → NDArray: Density of the first-passage time. Under a stress path it is the reference-stress density at the clock time tau(x) times the clock’s rate, AF of the stress in force at x.

ff(x: ArrayLike, Z: Any = None) → NDArray

Failure (CDF) of the first-passage time to the threshold.

For a model fitted with stress, Z is required: one stress row for a constant stress, or a StepSchedule for a stress profile; a model fitted with fit_from_df also takes a one-row DataFrame, read by column name. The same applies to every method below. A missing (nan) time or stress gives nan.

classmethod from_dict(model_dict: dict) → FirstPassageProcessModel: Rebuild a fitted process model from a to_dict() dict.

classmethod from_json(fp: str | PathLike) → Any: Load a model from a JSON file written by to_json(), or from the JSON text it returned (a string starting with {).

gamma: NDArray | None: Stress coefficients and the reference stress, for a model fitted with Z; both None otherwise.

hf(x: ArrayLike, Z: Any = None) → NDArray: Hazard function of the first-passage time.

property is_accelerated: bool: True for a model fitted with a stress Z.

mean(Z: Any = None) → float: Mean time to failure. At a constant stress it is the reference-stress mean divided by the acceleration factor; under a stress profile it is the integral of the survival function.

predict_rul(current_degradation: float, *, Z: Any = None, alpha_ci: float = 0.05) → ProcessRUL

Remaining useful life given the current degradation level.

Increments are independent for both processes, so the remaining first passage over the residual distance threshold - current_degradation follows the same law as a fresh process; its median and equal-tailed interval are returned.

Parameters:

current_degradation (float) – The unit’s current degradation level (not nan).
Z (array like or StepSchedule, optional) – For a model fitted with stress: the stress the unit will run at from now on – one row for a constant stress, or a StepSchedule whose time zero is now (or, for a model fitted with fit_from_df, a one-row DataFrame). It describes this one unit, so a missing value is refused.
alpha_ci (float, optional) – Tail probability of the returned interval, between 0 and 1. Default 0.05.

Returns:

The median remaining life and its equal-tailed interval.

Return type:

ProcessRUL

qf(p: ArrayLike, Z: Any = None) → NDArray: Quantile (inverse CDF) of the first-passage time (nan for a missing probability or stress).

random(size: int, Z: Any = None, random_state: int | None = None) → NDArray

Draw first-passage (failure) times from the fitted model.

Parameters:

size (int) – Number of draws.
Z (array like or StepSchedule, optional) – The stress, for a model fitted with Z (required then); each reference-stress draw is carried to calendar time along its clock.
random_state (int or numpy.random.Generator, optional) – Seed or generator for reproducible draws. None (the default) seeds from numpy’s global RNG, so np.random.seed controls it.

sf(x: ArrayLike, Z: Any = None) → NDArray: Survival function of the first-passage time.

to_dict() → dict: Serialise this fitted process model to a plain dict.

to_json(fp: str | PathLike | None = None, with_data: bool = False) → str | None

Write to_dict() to fp as strict JSON, or return it.

Parameters:

fp (str or os.PathLike, optional) – The file to write. Without it the JSON is returned as a string (as pandas.DataFrame.to_json does), which from_json also reads.
with_data (bool, optional) – Write to_dict(with_data=True), which also stores the fitted data, for the models whose to_dict takes with_data (the univariate Parametric and NonParametric); a TypeError for any other model. Defaults to False.

class surpyval.degradation.process_models.ProcessRUL(rul: float, rul_interval: tuple, prob_already_failed: float, alpha_ci: float)

Bases: object

Remaining-useful-life summary from a fitted process model.

rul

Median remaining useful life from the current state.

Type:: float

rul_interval

Equal-tailed 1 - alpha_ci interval for the remaining life.

Type:: tuple of float

prob_already_failed

1.0 when the current degradation is at or beyond the threshold (the remaining life and its interval are then 0), else 0.0.

Type:: float

alpha_ci

The tail probability of rul_interval.

Type:: float

Examples

predict_rul of a process model returns one. A unit that has degraded to 60 of a failure threshold of 100, drifting at 0.5 per unit time:

>>> from surpyval.degradation import WienerProcessModel
>>> model = WienerProcessModel(mu=0.5, sigma=1.0, threshold=100)
>>> rul = model.predict_rul(60.0)
>>> rul
ProcessRUL(rul=78.06, interval=(50.7, 120.3), prob_already_failed=0)
>>> round(rul.rul, 2)
78.06

Destructive Degradation

For tests that destroy the unit being measured, so each unit yields one observation at one time rather than a path. The degradation distribution at each time is modelled directly, and the failure distribution follows from the threshold crossing.

class surpyval.degradation.destructive.DestructiveDegradation_

Bases: object

Fitter for destructive degradation data (one destructive measurement per unit). Use the module-level singleton DestructiveDegradation.

fit(x: ArrayLike, y: ArrayLike, threshold: float, c: ArrayLike | None = None, distribution: ~typing.Any = <surpyval.univariate.parametric.distributions.lognormal.LogNormal_ object>, transform: str = 'linear', direction: str = 'auto') → DestructiveDegradationModel

Fit a destructive degradation model.

Parameters:

x (array_like) – Measurement time of each unit (one value per unit).
y (array_like) – The destructive degradation measurement of each unit.
threshold (float) – The degradation level D_f at which a unit is deemed failed.
c (array_like, optional) – Censoring of each measurement (not the time): 0 observed, 1 right-censored (e.g. did not break at the maximum load), -1 left-censored (below the test floor). Default all observed.
distribution (Parametric or str, optional) – Location-scale response distribution – LogNormal (default, positive response), Normal, or another such as Logistic or LogLogistic – as the object or its name. A distribution with positive support needs every measurement positive.
transform (str, optional) – Time transform \(\varphi(t)\) for the location: "linear", "log", "sqrt", "reciprocal", or "best" to pick the transform with the lowest AICc.
direction ({'auto', 'increasing', 'decreasing'}, optional) – Whether degradation moves up toward the threshold (wear) or down toward it (strength loss). "auto" infers it from the sign of the time-degradation trend.

Returns:

The fitted model, whose life-distribution methods (sf, ff, …) give the probability of having crossed threshold by each time.

Return type:

DestructiveDegradationModel

Warns:

UserWarning – “No finite maximum” when every measurement lies on the fitted path (noise-free readings): the fitted spread is then 0 to the precision of the fit, and sigma is meaningless.

Examples

Six units broken at each of four ages; strength falls log-linearly with age, and a unit has failed once its strength is below 20:

>>> import numpy as np
>>> from surpyval.degradation import DestructiveDegradation
>>> rng = np.random.default_rng(1)
>>> x = np.repeat([10.0, 20.0, 30.0, 40.0], 6)
>>> y = np.exp(4.0 - 0.02 * x + rng.normal(0, 0.1, 24))
>>> model = DestructiveDegradation.fit(x, y, threshold=20)
>>> model
Destructive Degradation Model
=============================
Response distribution : LogNormal
Time transform        : t
Direction             : decreasing
Threshold             : 20
Location              : 4.02814 + -0.0206616*t
Scale (sigma)         : 0.0612395
>>> model.sf([50, 80]).round(4)
array([0.4956, 0.    ])

fit_from_df(df: Any, x_col: str = 'x', y_col: str = 'y', c_col: str | None = None, **fit_kwargs: Any) → DestructiveDegradationModel

Fit a destructive degradation model from the columns of a pandas.DataFrame, with the argument names of DegradationAnalysis.fit_from_df.

Parameters:

df (pandas.DataFrame) – One row per unit tested.
x_col (str, optional) – Column of the measurement times. Defaults to "x".
y_col (str, optional) – Column of the measurements. Defaults to "y".
c_col (str, optional) – Column of the measurements’ censoring flags. Default all observed.
**fit_kwargs – Remaining arguments passed to fit(): threshold (required), and optionally distribution, transform and direction.

Returns:

The model fit() returns for the same arrays.

Return type:

DestructiveDegradationModel

Examples

>>> import numpy as np
>>> import pandas as pd
>>> from surpyval.degradation import DestructiveDegradation
>>> rng = np.random.default_rng(1)
>>> age = np.repeat([10.0, 20.0, 30.0, 40.0], 6)
>>> df = pd.DataFrame({
...     "age": age,
...     "strength": np.exp(4.0 - 0.02 * age + rng.normal(0, 0.1, 24)),
... })
>>> model = DestructiveDegradation.fit_from_df(
...     df, x_col="age", y_col="strength", threshold=20
... )
>>> model.sf([50, 80]).round(4)
array([0.4956, 0.    ])

class surpyval.degradation.destructive.DestructiveDegradationModel(distribution: Any, transform: str, direction: str, beta: NDArray, sigma: float, threshold: float, data: dict | None, neg_ll: float, transform_scores: dict | None = None)

Bases: SerialisableMixin

Result of DestructiveDegradation.fit().

Exposes the induced lifetime distribution at the failure threshold (sf / ff / Hf / df) plus the fitted degradation distribution over time (degradation_quantile). The fitted parameters are the location intercept and slope beta and the scale sigma.

Examples

Six units destroyed in a strength test at each of four ages; a unit has failed once its strength is below 20:

>>> import numpy as np
>>> from surpyval.degradation import DestructiveDegradation
>>> rng = np.random.default_rng(1)
>>> x = np.repeat([10.0, 20.0, 30.0, 40.0], 6)
>>> y = np.exp(4.0 - 0.02 * x + rng.normal(0, 0.1, 24))
>>> model = DestructiveDegradation.fit(x, y, threshold=20)

The median strength at ages 10 and 50, and the probability a unit is still above the threshold at 50 and 80:

>>> model.degradation_quantile(0.5, [10, 50]).round(3)
array([45.674, 19.987])
>>> model.sf([50, 80]).round(4)
array([0.4956, 0.    ])

Hf(x: ArrayLike) → NDArray: Cumulative hazard of the induced lifetime distribution.

cb(x: ArrayLike, on: str = 'sf', alpha_ci: float = 0.05, bound: str = 'two-sided', n_boot: int = 200, random_state: int | None = None) → NDArray

Bootstrap confidence bounds on the induced lifetime function on.

Units are resampled with replacement (each carrying its own (x, y, c)) and the whole fit is rerun, folding the estimation uncertainty into the band. Percentile bounds are returned.

Parameters:

x (array_like) – Times at which to evaluate the bound(s).
on ({'sf', 'ff', 'Hf'}, optional) – The lifetime function to bound ('R' and 'F' are accepted for 'sf' and 'ff'). Default 'sf'.
alpha_ci (float, optional) – Total tail probability. Default 0.05.
bound ({'two-sided', 'lower', 'upper'}, optional) – Two-sided bounds put [lower, upper] on the last axis, with alpha_ci / 2 in each tail. Default 'two-sided'.
n_boot (int, optional) – Number of bootstrap resamples. Default 200.
random_state (int or numpy.random.Generator, optional) – Seed or generator for the resampling. None (the default) seeds from numpy’s global RNG, so np.random.seed controls it.

degradation_quantile(p: ArrayLike, x: ArrayLike) → NDArray

The p-quantile of the destructive measurement at time x (the fitted degradation distribution dist(loc(x), sigma)).

Parameters:

p (float or array_like) – Probability (or probabilities) in (0, 1).
x (float or array_like) – Time(s) at which to read the degradation distribution; a scalar x gives a scalar result.

df(x: ArrayLike) → NDArray: Density of the induced lifetime distribution (finite-difference of the CDF; the closed form depends on the time transform).

ff(x: ArrayLike) → NDArray: Failure (CDF) of the lifetime induced by crossing the threshold.

classmethod from_dict(d: dict) → DestructiveDegradationModel

Rebuild a model from a to_dict() dictionary.

Dictionaries written before the fit data was stored still load; the model they give predicts, but its cb() raises because there is no data to resample.

See also

to_dict, to_json, from_json

classmethod from_json(fp: str | PathLike) → Any: Load a model from a JSON file written by to_json(), or from the JSON text it returned (a string starting with {).

median_degradation(x: ArrayLike) → NDArray: Median destructive measurement at time x.

sf(x: ArrayLike) → NDArray: Reliability of the induced lifetime distribution.

to_dict() → dict

Serialise this fitted model to a plain, JSON-serialisable dict.

The fit data (x, y, c) is stored along with the fitted parameters – as DegradationModel stores its raw data – so the restored model reproduces the original’s predictions and its bootstrap cb() (with the same random_state, exactly).

See also

from_dict, to_json, from_json

to_json(fp: str | PathLike | None = None, with_data: bool = False) → str | None

Write to_dict() to fp as strict JSON, or return it.

Parameters:

fp (str or os.PathLike, optional) – The file to write. Without it the JSON is returned as a string (as pandas.DataFrame.to_json does), which from_json also reads.
with_data (bool, optional) – Write to_dict(with_data=True), which also stores the fitted data, for the models whose to_dict takes with_data (the univariate Parametric and NonParametric); a TypeError for any other model. Defaults to False.