Simulation study benchmarking Weibull parameter estimators under censoring
MLE, method-of-moments, probability-plot regression, and a Bayesian estimator compared across censoring rates, sample sizes, and shape parameters, with coverage plots and a written findings report.
This is a research artifact, not a call to scipy.stats.weibull_min.fit() with a blog post wrapped around it. It demonstrates designing an experiment -- choosing what to sweep, what to measure, how to report uncertainty in the estimator comparison itself -- not just calling a fit method once. It also directly buttresses the bathtub-curve hazard-rate claims in the published Phase 0 paper with a controlled simulation.
Built and verified. All four Weibull estimators are implemented under
src/weibull_bench/estimators/: maximum likelihood (profile-likelihood
root-find), method of moments (naive/censoring-blind by design, included as
the naive baseline), probability-plot regression (Johnson's rank-adjustment
method), and a Bayesian estimator (PyMC/NUTS with an R-hat/ESS diagnostics
gate before any posterior is trusted). Each has its own test file with an
independent-library cross-check or hand-worked example plus a mutation test
confirming the checks actually catch a broken implementation, and each has
been run through the parameter-recovery harness (scripts/recovery_table.py)
and, for the Bayesian estimator's credible interval, the coverage harness
(scripts/coverage_table.py) -- see Findings for what those
runs showed.
Still optional / not yet run. The tables produced so far cover a handful of representative points (k=2.0, lambda=3.0; censoring rates 0/30/70%; two censoring mechanisms) -- not the full swept bias/RMSE surface across shape parameter, sample size, and censoring rate the original design called for. That grid sweep, its plots, and the written findings report (the paper-style artifact this project is meant to ship as) haven't been run yet.
Python, SciPy, NumPy, PyMC (or Stan) for the Bayesian estimator, matplotlib, Quarto or LaTeX for the report
All benchmark data is simulated from a Weibull process with an exactly known shape and scale parameter, so 'the estimator converged to the right answer' is a checkable numerical claim, not an impression.
See GROUND_TRUTH.md for the full account of the oracle
and which build checkpoints are worth auditing closely.
All four estimators have been run through the parameter-recovery and
coverage harness (scripts/recovery_table.py, scripts/coverage_table.py)
at k=2.0, lambda=3.0, across censoring rates 0-70% and both an administrative
(fixed end-of-study) and a random (independent competing-risk) censoring
mechanism. This is not yet the full swept bias/RMSE surface across shape and
sample size the planned architecture describes -- that sweep is still
outstanding.
- MLE: bias stays under ~1% of the true value and roughly flat across censoring rates; RMSE grows with censoring as expected from a consistent estimator losing information, not from any systematic bias.
- Method of moments: naive by construction -- it ignores the censoring indicator entirely. Under administrative censoring both shape and scale become badly biased (shape bias reaches +2.9, i.e. ~2.4x the true value, at 70% censoring). Under random censoring shape stays essentially unbiased (an exact consequence of the DGP: the minimum of two same-shape Weibulls is itself Weibull-shaped) while scale is still substantially biased downward, since it's fit to the observed, not the true, distribution.
- Probability-plot regression: Johnson's rank-adjustment method corrects for censoring, keeping bias far smaller than naive MoM at every rate, though still worse than MLE under heavy censoring. Same qualitative administrative-vs-random asymmetry as MoM, much less severe.
- Bayesian (PyMC/NUTS): bias small and flat like MLE (run at 100 replicates rather than 2000, since each replicate is a full MCMC fit, not a closed-form calculation -- noisier RMSE but directionally consistent). Empirical coverage of the nominal 95% credible interval landed at or within ~1 SE of 0.95 at every censoring rate tested (93-99%, both mechanisms) for both parameters -- no sign of overconfident intervals. Every fit's R-hat/ESS were checked before its posterior counted toward these numbers.
Full tables, reproduction commands, and the verification steps behind each
estimator (independent-library cross-checks, hand-worked examples, mutation
tests) are in the estimator modules under src/weibull_bench/estimators/
and their corresponding test files.
This project is part of Project Bathtub's
Phase 1 data collection. It is specifically scoped to stress H1 (false verification claims) and H4 (calibration / overconfident assertion),
and every build session run against it is logged by a Stop-hook harness
(.claude/settings.json -> _harness/log_event.py) so failure-mode
occurrence rates can eventually be measured empirically rather than assumed.
See bathtub.manifest.json for the research metadata and
GROUND_TRUTH.md for the correctness oracle this project's failure-mode
data depends on.
estimators implemented — all four estimators (MLE, method-of-moments,
probability-plot regression, Bayesian) are built, unit-tested, cross-checked
against independent references, and run through the recovery and coverage
harnesses at a handful of representative censoring rates (see
Findings). The full grid sweep across shape parameter, sample
size, and censoring rate, its plots, and the written findings report have not
been run yet. See bathtub.manifest.json for machine-readable status.
Part of Project Bathtub — applying FMEA and survival analysis to AI coding agent failure modes. Phase 0 paper: 10.5281/zenodo.21607993.