I have been working on the same statistical problem since 2020. I am now a PhD student in CS. The problem has not changed, but my understanding of it has, and the tools I have built around it look nothing like what I started with.

The problem: a series system fails when any component fails. You observe system-level failure times. But you often cannot tell which component caused the failure (masking). Some systems are still running when testing ends (censoring). Given this incomplete data, estimate component reliability.

This is not a tutorial. It is a map of where things stand and where they are going.

Last week I announced maskedcauses, the R package for estimating component reliability from masked series system failures. That post covered the three likelihood models and the path to CRAN.

This post is about what happened next: the package now supports four observation types (exact, right-censored, left-censored, and interval-censored) via composable observation functors. Along the way, I wrote four vignettes, removed the md.tools dependency, and developed a verification methodology for keeping prose honest about simulation results.

Note (February 2026): This package has been renamed from likelihood.model.series.md to maskedcauses.

Two days ago, I submitted likelihood.model to CRAN, the foundation package for composable statistical inference. Next in line: maskedcauses, which implements maximum likelihood estimation for series systems where component failure causes are masked.

This package is the practical result of my master’s thesis work. Three years of theoretical development, now packaged for anyone analyzing masked failure data.

The Problem: Masked Component Failures

A series system fails when any of its \(m\) components fails. In reliability testing, you observe the system fail at time \(t\), but two layers of uncertainty obscure the full picture:

1. Right-censoring: Some systems are still running when testing ends. You know they survived at least until time \(\tau\), but not how much longer they would have lasted.

2. Masked cause of failure: When a system fails, you often can’t identify which component caused it. Diagnostic tests might narrow it down to a candidate set of possible causes, but the true failure component remains ambiguous.

This happens constantly in practice. Electronic systems fail with only board-level diagnostics. Industrial machinery fails without root-cause teardown. Medical devices fail with symptoms pointing to multiple possible subsystems.

The question: given this incomplete information, can you still estimate the lifetime distribution of each component?

The Package: Three Likelihood Models

maskedcauses provides three models with different complexity-accuracy tradeoffs:

Each model implements the full inference stack: loglik(), score(), hess_loglik(), rdata(), and assumptions().

The C1-C2-C3 Conditions

The models assume three conditions that simplify the likelihood:

• C1: The failed component is in the candidate set with probability 1
• C2: Given the failed component is in the candidate set, masking probability is uniform across candidates
• C3: Masking probabilities are independent of system parameters \(\theta\)

Under these conditions, the masking mechanism factors out of the likelihood. You can estimate component parameters without modeling the diagnostic process itself. That’s why the package name includes “c1_c2_c3”.

symlik is a Python library for symbolic likelihood models. Write your log-likelihood as a symbolic expression, and it derives everything needed for inference.

The Problem

Traditional statistical computing gives you two choices:

1. Manual derivation. Work out score functions and information matrices by hand, then implement them. Error-prone, tedious.
2. Numerical approximation. Use finite differences. Unstable, slow, no symbolic form to inspect.

The Approach

symlik takes a third path: symbolic differentiation. Define the model once, get exact derivatives automatically.

from symlik.distributions import exponential

model = exponential()
data = {'x': [1.2, 0.8, 2.1, 1.5]}

mle, _ = model.mle(data=data, init={'lambda': 1.0})
se = model.se(mle, data)

print(f"Rate: {mle['lambda']:.3f} +/- {se['lambda']:.3f}")
# Rate: 0.714 +/- 0.357

Behind the scenes, symlik:

1. Symbolically differentiates the log-likelihood to get the score function
2. Differentiates again for the Hessian
3. Computes Fisher information from the Hessian
4. Derives standard errors from the inverse information matrix

All exact. No numerical approximation.

Custom Models

The real power is defining custom models using s-expressions:

from symlik import LikelihoodModel

# Exponential: l(lambda) = sum[log(lambda) - lambda*x_i]
log_lik = ['sum', 'i', ['len', 'x'],
           ['+', ['log', 'lambda'],
            ['*', -1, ['*', 'lambda', ['@', 'x', 'i']]]]]

model = LikelihoodModel(log_lik, params=['lambda'])

# Symbolic derivatives available
score = model.score()       # Gradient
hess = model.hessian()      # Hessian matrix
info = model.information()  # Fisher information

You define the log-likelihood once as a symbolic expression. symlik computes the rest.

Heterogeneous Data

One of symlik’s strengths is handling mixed observation types, which is exactly what you need for reliability analysis with censored data:

from symlik import ContributionModel
from symlik.contributions import complete_exponential, right_censored_exponential

model = ContributionModel(
    params=["lambda"],
    type_column="status",
    contributions={
        "observed": complete_exponential(),
        "censored": right_censored_exponential(),
    }
)

data = {
    "status": ["observed", "censored", "observed", "observed", "censored"],
    "t": [1.2, 3.0, 0.8, 2.1, 4.5],
}

Each observation type contributes differently to the likelihood. symlik handles the bookkeeping.

Connection to Research

symlik is the Python successor to my R package likelihood.model. It implements the theoretical framework from my thesis work on likelihood-based inference for series systems.

The Weibull Series Model Selection paper shows applications to reliability engineering, the kind of complex likelihood that benefits from symbolic treatment.

Powered by rerum

symlik uses rerum for symbolic differentiation. rerum is a pattern matching and term rewriting library that handles the calculus. The separation means you can use rerum for other symbolic computation tasks beyond likelihood models.

Installation

Available on PyPI:

pip install symlik

Documentation at queelius.github.io/symlik.

See the project page for more details.

In a series system, the system fails when any component fails. You observe the system failure time \(t\) and a candidate set \(C \subseteq \lbrace 1,2,\ldots,m\rbrace\) of components that might have caused the failure. But you do not know which component in \(C\) actually failed. This is masked failure data.

The standard approach is numerical optimization of the likelihood. This paper shows that for exponential component lifetimes, everything has a closed form.

Closed-Form Fisher Information

For exponential masked data with arbitrary masking patterns:

where \(\hat{\omega}_A\) is the observed frequency of candidate set \(A\). You can compute asymptotic variances directly, check identifiability before running any estimation, and analyze optimization stability. All without fitting a model first.

Sufficient Statistics

The mean system lifetime and the candidate set frequency vector are sufficient statistics. That reduces an entire dataset to \(1 + \binom{m}{w}\) numbers, where \(w\) is the masking width.

This is a real simplification. All the statistical information in your data is captured by two things: how often each candidate set appears, and what the average failure time is. Nothing else matters for inference.

Closed-Form MLE for Three Components

For \(m=3\) components with pairwise masking (\(w=2\)), the MLE has an explicit closed-form solution:

No numerical optimization. No iterative algorithms. Just plug in your sufficient statistics.

The \(w=2\) case is the interesting one. \(w=1\) means no masking (you know exactly which component failed). \(w=m\) means complete masking (the candidate set is always everything, so you have no diagnostic information). \(w=2\) is the simplest case where masking actually matters, and it is the one where closed-form solutions exist.

Asymptotic Theory

The MLE follows:

with explicit Wald-type confidence intervals using the closed-form Fisher information. So you get point estimates and uncertainty quantification, all analytically.

Why Exponential?

The exponential assumption is not just for tractability, though it helps. Constant hazard rate models systems subject to random external shocks. The memoryless property simplifies the likelihood structure. And exponential is the foundation for generalization to Weibull and other distributions.

Most R packages hardcode specific likelihood models. likelihood.model takes a different approach. Likelihoods are first-class objects that compose, and the framework is generic enough to work with any distribution.

The Interface

A likelihood model is anything implementing these generic methods:

• loglik(model, data, params) – log-likelihood
• score(model, data, params) – score function (gradient)
• hessian(model, data, params) – observed information matrix

That is the interface. If your model implements these three methods, it plugs into the entire MLE stack: optimization, confidence intervals, hypothesis testing, model selection. You do not couple to specific distributions.

Likelihood Contributions

The key class is likelihood_contr_model, a likelihood built from independent contributions:

# Different observation types get different likelihood contributions
model <- likelihood_contr_model(
  exact = normal_contrib(),
  right_censored = censored_contrib()
)

This handles heterogeneous data in a unified framework. You can mix exact observations, right-censored observations, truncated observations, and different distribution families within one model. Each observation type gets its own likelihood contribution, and they combine additively in log-space.

Why This Design

The i.i.d. assumption decomposes a joint likelihood into additive log-likelihood contributions. That is how MLE actually works. likelihood.model makes this decomposition explicit and compositional.

Likelihood models are objects you manipulate, not function calls buried inside a fitting routine. You can build complex models from simple, independent pieces. You can swap in different contribution types without rewriting the rest of your code. And because the interface is generic, it works with algebraic.mle for fitting, hypothesize for testing, and any optimization backend that speaks the same protocol.

This is the same compositional philosophy as my thesis work on masked failure data. Series systems with masked causes have multiple observation types (masked vs. unmasked, different candidate sets) that each contribute differently to the likelihood. likelihood.model handles that naturally.

R package – MIT licensed – Documentation – GitHub

Maximum likelihood estimators have rich mathematical structure. They are consistent, asymptotically normal, efficient. algebraic.mle exposes this structure through an algebra where MLEs are objects you compose, transform, and query.

The Abstraction

An MLE is not just a vector of parameter estimates. It is a statistical object that carries point estimates \(\hat{\theta}\), the Fisher information matrix \(I(\hat{\theta})\), the variance-covariance matrix \(I^{-1}(\hat{\theta})\), Wald-type confidence intervals from asymptotic normality, the log-likelihood value, and convergence diagnostics.

The package wraps all of this in a consistent interface:

library(algebraic.mle)

fit <- mle(likelihood_model, data)
coef(fit)           # Parameter estimates
vcov(fit)           # Variance-covariance matrix
confint(fit)        # Confidence intervals
logLik(fit)         # Log-likelihood
aic(fit)            # Model selection

Composition

The real point is that MLEs compose. Independent models combine:

fit1 <- mle(model1, data1)
fit2 <- mle(model2, data2)
combined <- fit1 + fit2  # Joint likelihood

The package handles the algebra. Joint log-likelihood, block-diagonal Fisher information, everything propagates correctly. This works because likelihoods from independent data sources multiply, and multiplication of likelihoods is addition of log-likelihoods. That is a monoid. The package enforces it.

The Ecosystem

algebraic.mle is the foundation for a family of packages:

The typical workflow:

1. Define distributions with algebraic.dist
2. Specify likelihood contributions with likelihood.model
3. Fit the model and get an mle object from algebraic.mle
4. Query statistical properties: confidence intervals, hypothesis tests, model selection

For series systems with masked data:

library(maskedcauses)
library(algebraic.mle)

# Specify masking model (C1-C2-C3 conditions)
model <- md_likelihood_model(components = 3, masking = "bernoulli")

# Fit -> returns algebraic.mle object
fit <- md_mle_exp_series_C1_C2_C3(masked_data)

# All the standard MLE methods work
confint(fit)
vcov(fit)
aic(fit)

Theory

The asymptotic properties that algebraic.mle exploits come from classical MLE theory:

The expo-masked-fim paper derives closed-form Fisher information for exponential series systems. That is exactly what algebraic.mle uses internally for variance estimation in that case.

For more complex models (Weibull, relaxed masking conditions), we compute Fisher information numerically via observed information:

Design Principles

Separation of concerns. The likelihood specification (likelihood.model) is independent of the fitting algorithm (numerical.mle) and the result type (algebraic.mle). You can swap optimizers without changing downstream code.