Parametric Recurrence Model

class surpyval.recurrent.parametric.parametric_recurrence.ParametricRecurrenceModel

Bases: RecurrenceSimulationMixin, LikelihoodInferenceMixin

A class for holding the parameters, data, and usefult methods for a fitted parametric recurrence model. This is the result of the fit calls from the counting distributions.

When fitted by maximum likelihood the model also carries the likelihood- inference behaviour (log_likelihood, aic, bic, standard_errors) from LikelihoodInferenceMixin. Models built by from_params or fitted by how="MSE" carry no likelihood, so those methods raise.

Example

>>> from surpyval import Exponential
>>> from surpyval.recurrent import HPP
>>> import numpy as np
>>> np.random.seed(1)
>>> x = Exponential.random(10, 1e-3).cumsum()
>>> model = HPP.fit(x)

cif(x)

Compute the cumulative incidence function (CIF) based on the fitted model. No need to pass parameters as it uses the parameters of the fitted model.

Parameters: x (array_like) – Values at which to compute the CIF.
Returns: Computed cumulative intensity function values.
Return type: array_like

cif_cb(x, alpha_ci=0.05, bound='two-sided')

Confidence bounds on the fitted CIF at x, from the delta method.

The variance of the fitted CIF is propagated from the parameter covariance (the inverse observed information) through the CIF’s gradient, and the bounds are computed on the log scale – the same construction as the exponential Greenwood bounds on the nonparametric MCF – so they cannot go negative.

Parameters

x (array_like) – Values at which to compute the confidence bounds.
alpha_ci (float, optional) – The total tail probability of the bound(s). Default is 0.05.
bound ({'two-sided', 'lower', 'upper'}, optional) – Two-sided bounds are returned as an (len(x), 2) array with columns [lower, upper]; one-sided bounds have the shape of x.

Returns

The confidence bounds on the CIF.

Return type

numpy array

count_terminated_simulation(events, items=1, seed=None)

Simulate count-terminated recurrence data based on the fitted model.

Parameters

events (int) – Number of events to simulate.
items (int, optional) – Number of items (or sequences) to simulate. Default is 1.
seed (int or numpy.random.Generator, optional) – Seed for a reproducible simulation. When None (default) the numpy global RNG is used.

Returns

An NonParametricCounting model built from the simulated data.

Return type

NonParametricCounting

count_terminated_simulation_data(events, items=1, seed=None)

Simulate count-terminated recurrence data and return the raw events.

Unlike count_terminated_simulation() (which returns the fitted NonParametricCounting MCF), this returns the simulated event data itself, ready to be refitted or inspected via .x/.i/.c/ .n.

Parameters

events (int) – Number of events to simulate per sequence.
items (int, optional) – Number of items (or sequences) to simulate. Default is 1.
seed (int or numpy.random.Generator, optional) – Seed for a reproducible simulation.

Returns

The simulated recurrence data in xicn format.

Return type

RecurrentEventData

Notes

Count termination is a failure-terminated (Type II) scheme: each item is observed until its events + 1-th event, so its observation window is the random time of that last event and every event is exact (c = 0). Parametric fits handle this correctly – the interarrival/intensity likelihood ends at the last observed event and the MLE is consistent. The nonparametric MCF, however, is only reliable up to roughly events recurrences: beyond that the at-risk set is depleted and the curve is biased (which is why count_terminated_simulation() trims to mcf_hat < events). For a fixed-window observation scheme, use time_terminated_simulation_data(), which right-censors each item at T.

covariance(): Approximate parameter covariance matrix, ordered to match parameter_names. Computed as the inverse of the numerical Hessian of the negative log-likelihood at the MLE.

cramer_von_mises(n_boot=200, seed=None)

Cramer-von Mises goodness-of-fit test of the fitted intensity.

Conditional on the number of events an item shows in its observation window, the transformed times [cif(t) - cif(entry)] / [cif(close) - cif(entry)] are iid U(0, 1) when the fitted intensity is the true one; the Cramer-von Mises statistic measures their departure from uniformity (for the power-law process this is the construction behind Crow’s goodness-of-fit test). Because the parameters were estimated from the same data, the p-value is computed by a parametric bootstrap: data is simulated from the fitted model over the same observation windows, refitted, and the statistic recomputed. Failure-truncated items (no explicit window close) are approximated in the bootstrap by a window fixed at their last observed event.

Parameters

n_boot (int, optional) – Number of bootstrap replicates for the p-value. Default is 200.
seed (int or numpy.random.Generator, optional) – Seed for a reproducible p-value.

Returns

The observed statistic and its bootstrap p-value.

Return type

GoodnessOfFitResult

iif(x)

Compute the intensity function based on the fitted model. No need to pass parameters as it uses the parameters of the fitted model.

Parameters: x (array_like) – Values at which to compute the intensity.
Returns: Computed instantaneous intensity functions values.
Return type: array_like

mcf(x)

The mean cumulative function (MCF). For these counting processes the MCF equals the cumulative intensity, so this is a closed-form alias for cif() (overriding the simulation-based estimate in the mixin).

Parameters: x (array_like) – Values at which to compute the MCF.
Returns: The MCF evaluated at x.
Return type: array_like

param_cb(name, alpha_ci=0.05, bound='two-sided')

Confidence bound(s) on a fitted parameter, mirroring the univariate Parametric.param_cb API.

Wald bounds from the observed information, computed on a transformed scale chosen from the parameter’s bounds so the result respects its support: log scale for one-sided-bounded parameters (e.g. a positive rate), logit scale for interval-bounded parameters (e.g. a repair efficiency in (0, 1)), and the natural scale for unbounded ones.

Parameters

name (str) – The parameter to bound; one of parameter_names.
alpha_ci (float, optional) – The total tail probability of the bound(s). Default is 0.05.
bound ({'two-sided', 'lower', 'upper'}, optional) – Two-sided bounds are returned as [lower, upper].

Returns

The confidence bound(s) on the parameter.

Return type

numpy array

plot(ax=None, plot_bounds=True, confidence=0.95)

Plot the fitted CIF over the nonparametric MCF of the data used to fit it, with a delta-method confidence band around the fitted curve when the model carries a likelihood.

Parameters

ax (matplotlib axes, optional) – An axes object to draw the plot on. Creates a new one if not provided.
plot_bounds (bool, optional) – Whether to draw the confidence band around the fitted CIF. Ignored for models with no likelihood (how="MSE" fits and from_params models). Default is True.
confidence (float, optional) – The confidence level of the band. Default is 0.95.

Returns

An axes object with the plot.

Return type

matplotlib axes

residuals(kind='cumulative_hazard')

Residual diagnostics for the fitted model, from the time-rescaling theorem.

Parameters: kind ({'cumulative_hazard', 'pit', 'martingale'}, optional) – 'cumulative_hazard' returns the rescaled interarrival times cif(t_k) - cif(t_{k-1}) of every observed event (pooled across items), which are iid Exp(1) under the fitted model. 'pit' applies the probability integral transform 1 - exp(-e) to those residuals, giving iid U(0, 1) values. 'martingale' returns one residual per (sorted-unique) item: its observed event count minus the count the model expects over its observation window; positive values mean the item saw more events than predicted.
Returns: The residuals.
Return type: numpy array

standard_errors(): Standard errors of the fitted parameters (the square roots of the diagonal of covariance()), ordered to match parameter_names. Entries are NaN where the variance is non-positive, which typically indicates a boundary optimum.

time_terminated_simulation(T, items=1, tol=1e-08, max_events=10000, seed=None)

Simulate time-terminated recurrence data based on the fitted model.

Parameters

T (float) – Time termination value.
items (int, optional) – Number of items (or sequences) to simulate. Default is 1.
tol (float, optional) – Interarrival times below this value end the sequence early; a tiny increment indicates the cumulative time has stalled below T (a possible asymptote). Default is 1e-8.
max_events (int, optional) – Hard cap on the number of events simulated per sequence. This is the backstop that guarantees termination for sequences whose cumulative time cannot reach T. Default is 10000.
seed (int or numpy.random.Generator, optional) – Seed for a reproducible simulation. When None (default) the numpy global RNG is used.

Returns

An NonParametricCounting model built from the simulated data.

Return type

NonParametricCounting

Warning

A sequence is terminated early and right-censored at its last event if an interarrival time falls below tol or it reaches max_events before T. A warning is raised in either case.

time_terminated_simulation_data(T, items=1, tol=1e-08, max_events=10000, seed=None)

Simulate time-terminated recurrence data and return the raw events.

Unlike time_terminated_simulation() (which returns the fitted NonParametricCounting MCF), this returns the simulated event data itself, ready to be refitted or inspected via .x/.i/.c/ .n. Each sequence is right-censored at T.

Parameters

T (float) – Time termination value.
items (int, optional) – Number of items (or sequences) to simulate. Default is 1.
tol (float, optional) – Interarrival times below this value end the sequence early. Default is 1e-8.
max_events (int, optional) – Hard per-sequence event cap that guarantees termination. Default is 10000.
seed (int or numpy.random.Generator, optional) – Seed for a reproducible simulation.

Returns

The simulated recurrence data in xicn format.

Return type

RecurrentEventData

trend_test(test='laplace', alternative='two-sided')

Run a trend test on the data this model was fitted to.

This is a convenience wrapper around the standalone tests in surpyval.recurrent.tests – the null hypothesis is that the events follow a homogeneous Poisson process (no trend), so it checks whether the data warranted a time-varying intensity at all. The model’s parameters play no part in the statistic.

Parameters

test ({'laplace', 'mil_hdbk_189c'}, optional) – The trend test to run. Default is ‘laplace’.
alternative ({'two-sided', 'increasing', 'decreasing'}, optional) – The alternative hypothesis. Default is ‘two-sided’.

Returns

The test result, carrying the statistic, p-value and suggested trend direction.

Return type

TrendTestResult