Changelog

All notable changes to this project are documented in this file.

The format is based on Keep a Changelog, and this project adheres to Semantic Versioning.

Unreleased

0.6.3 - 2026-08-30

Cut so that the corrected replication package reaches an archived record. 0.6.2 archived a 31-file package; #409 and #411 changed it materially, so the DOI the paper cited no longer described the package in the repository.

Fixed

• pyproject.toml claimed the released version while main carried commits past the tag, so a wheel built from main would have announced itself as an already published release from a different tree. After a release main now moves to the next .devN while CITATION.cff goes on naming the last release.
• build_online_resource.py printed a note and exited 0 when the DOI or the corresponding e-mail was still a placeholder, after writing a complete-looking journal supplement. It refuses before writing anything now, with --allow-placeholders for drafts.

Changed

• The replication notes state two asymmetries rather than leaving them to be discovered: the frozen stress artifact carries no software provenance, and the stress layer ships no replicate-level export, so its Monte Carlo standard errors can be read but not independently recomputed.
• ROADMAP.md and README.md describe the package that exists. Both listed as future work things that ship, and the README's module table named five modules of thirteen.

0.6.2 - 2026-08-30

Added

• The replication package for the sparse-contamination paper is now tracked at research/sparse_contamination/replication_package/: simulation drivers, analysis-only scripts, the frozen summary and per-replicate losses, the configuration and provenance records, the manuscript with its generated fragments, and a SHA-256 manifest over all 31 files. The paper promised a deposited archive; keeping it here means the release tarball Zenodo archives contains it, so it is citable without a separate deposit. It is excluded from the PyPI distributions, which ship only the package, tests and two metadata files.

• make verify-release VERSION=x.y.z, backed by scripts/verify_release.py, asks whether a released version is actually public: tag, GitHub release, PyPI and the Zenodo archive, plus the local metadata that names it. The preflight cannot answer this, because half of a release does not exist when it runs -- Zenodo mints the version DOI only on archiving. Run against 0.6.0 it reports the incident exactly: tagged, released on GitHub, absent from PyPI.

• make release-check and make release-preflight, backed by scripts/check_release.py, verify that the five files carrying the release identity agree before a release is tagged. The version lives in pyproject.toml, CITATION.cff, CHANGELOG.md, docs/about/citation.md and README.md; the date lives in CITATION.cff, .zenodo.json and CHANGELOG.md. Most of those couplings already had tests, but they ran in CI, and the publish job depends on the test job -- so a release commit that missed a file still tagged cleanly and simply never reached PyPI. The changelog was checked by nothing at all. --pre-tag additionally refuses a version whose tag is already spent, since PyPI will not accept a re-upload.

• docs/development/releasing.md now lists all five files. It previously named only pyproject.toml and CITATION.cff, which is why 0.6.0 was released with a stale changelog and stale citation guidance.

Fixed

• The APA, IEEE, MLA and Chicago citations in docs/about/citation.md, and the suggested citation in README.md, each said "Version 0.6.1" beside the concept DOI, which resolves to whatever is newest. Five citations named a version they did not pin. check_release.py now reads citations by block -- a fenced snippet or a run of blockquote lines, since the version and the DOI are rarely on the same line -- and requires any block naming a version to cite that version's DOI.
• Two false statements in docs/development/releasing.md. coverage does not run "on releases and on main"; it runs on a release or a manual dispatch and on nothing else. And the After Release checklist still said publish depends on test alone, which the Before Tagging section had already been corrected away from in the same file.

Changed

• Releasing now dispatches coverage before tagging. It is the one gate no pull request exercises, publish waits on it, and on an ordinary release it first runs against the shipped code after the tag is public -- so a failure there costs a version number for the same reason 0.6.0 did.

• The "Citing release X exactly" block in docs/about/citation.md quoted the concept DOI, which by that same page's table "always resolves to the most recent release". The block that exists to pin a version pinned nothing. It now carries the version DOI, and check_release.py fails when the DOI named in CITATION.cff is not the one the citation guidance quotes. The gap survived because a version DOI does not exist when a release is cut -- Zenodo mints it on archiving, and updating the docs afterwards was a step nothing enforced.

Changed

• CITATION.cff now names the version DOI for the release it describes, and lists all seven minted DOIs. The three that were pending have been minted: 0.6.1 is 10.5281/zenodo.22166257, 0.6.0 is 10.5281/zenodo.22161110 and 0.5.0 is 10.5281/zenodo.22070996.
• docs/development/releasing.md records that publish waits on five jobs, not one, and that coverage runs only on releases and main -- so no pull request exercises that gate before a release depends on it.

0.6.1 - 2026-08-29

Fixed

• .zenodo.json carried publication_date: 2026-08-23 while CITATION.cff had been moved to 2026-08-29. scripts/validate_zenodo_metadata.py requires the two to agree, so the CI test job failed on every Python version and the release job, which depends on it, never published 0.6.0 to PyPI. The dates now agree. There is no code change: 0.6.1 is 0.6.0 with metadata that passes its own validator.

0.6.0 - 2026-08-29

Fixed

• The discrete distribution functions were not zero below their support. Zipf.cdf clamped its argument up to 1 before summing, so cdf(0) -- and cdf(-5) -- returned P(X = 1), which for Zipf(s=1.5) is 0.387 rather than 0. DiscretePareto.cdf did the same at k_min: with k_min = 3, cdf(2) returned 0.165. Clamping the upper end is right, since the support is truncated there; clamping the lower end is not, because below the support the answer is zero and not the first atom.

Changed

• Zipf, YuleSimon and DiscretePareto take a number or an array and mirror what they were given, as the continuous families have since 0.5.0. ppf returns an integer, or an integer array.
• Zipf.ppf and DiscretePareto.ppf find the quantile by searching a cumulative table instead of walking the support one k at a time. The scan was fine near the median and ruinous in the tail, which is where this library is used: at Zipf(s=1.1, kmax=100_000), 200 calls to ppf(0.99) took 2.31s and now take 2.0ms; 20,000 quantiles took 9.24s and now take 2.0ms. Values are unchanged -- the same k for every probability tested, and seeded samples identical.

Added

• heavytails.registry: one place mapping a name to a distribution (#383). Names resolved in at least four ways before — a table in cli.py, getattr(heavytails, name) in performance.py and roadmap.py, and if/elif chains in validation.py and utilities.py, the chain repeated five times in validation.py alone. They disagreed: the CLI took student-t and gpd, the chains took studentt and had no entry for the generalized Pareto, getattr took only the class name. Every name any of them accepted is a registry alias, so nothing that worked stops working, and the CLI's table is now a view of the registry rather than a second list.

Changed

• heavytails.performance, extensions, utilities and validation are type-checked (#384). Their ignore_errors suppressions are gone. Most of the 47 errors were the lookup chains rebinding one variable to five distribution types; the registry removed those, and the rest were missing annotations.
• PerformanceBenchmarks.benchmark_sampling declares dict[int, dict[str, Any]]. It returns one entry per sample size, keyed by that size, and previously claimed dict[str, Any].

Removed

The 59 shipped functions that raised NotImplementedError (#312). A module of empty signatures is indistinguishable from a module of working code until you call it, and these set an expectation the package does not meet.

• heavytails.extensions: RegimeSwitchingModel, VineCopula, SpatialHeavyTailProcess, ScientificComputingIntegration, BayesianHeavyTailRegression, DistributionClassifier, DomainSpecificExtensions, TailGeometryAnalysis, and HeavyTailCopula.sample. The working copulas there are untouched.
• heavytails.utilities: ConfigurationManager, DataQualityAssessment, DistributionMetadata, FinancialDataScraper, CitationManager, InteractiveTutorials, PluginManager, UnitConverter. DataIO, AutoFit, ParameterValidator and StatisticalSummary are untouched.
• heavytails.roadmap: MultivariateStudentT, HeavyTailGARCH, HeavyTailSurvival, HeavyTailsAPI, gpu_sampling, adaptive_threshold_selection, bayesian_parameter_estimation, memory_efficient_sampling, diagnostic_plots. The stub MultivariateStudentT shadowed the working one in heavytails.multivariate, which is the clearest case of the harm.
• heavytails.performance: cython_special_functions, jit_accelerated_functions, optimized_rejection_sampling, distribution_specific_optimizations, MemoryProfiler, OnlineEstimation.
• heavytails.validation: fuzz_testing, special_function_accuracy_analysis, python_version_compatibility, ParameterEstimationValidation, RegressionTesting, MathematicalPropertyVerification.

The ideas are recorded in ROADMAP.md under "Considered, Not Committed", where they promise nothing. The distribution registry is kept as work rather than dropped.

Eighteen tests went with them, each of which asserted only that a stub raised NotImplementedError.

Fixed

• heavytails.viz.plot_tail overlaid a sample drawn from the fitted model rather than the model's survival curve. That put Monte Carlo noise into the reference the reader compares against, worst in the far tail where its last points rested on a handful of draws: against Pareto(alpha=2) at n=1000 the reference wandered up to 1.238 in log survival from the curve it claimed to be -- a factor of three, seed-dependent, and entirely in the region the plot exists to show. It reads as misfit and is not. The overlay is now fitted.sf evaluated on a grid, which is exact, monotone and independent of any seed. Points where a bounded model has zero survival are not drawn, since log(0) is not a point on the plot.

0.5.0 - 2026-08-23

NumPy is now required, and every distribution method takes a number or an array. The library was pure Python by choice, and that choice stopped paying: evaluating a density over a hundred thousand points cost a hundred thousand interpreter round trips. Over a million points it now takes about 20 milliseconds.

That is the reason this is 0.5.0 and not 0.4.1. Installing into an environment that cannot have NumPy will no longer work.

Breaking

• NumPy is a runtime dependency (numpy>=1.24,<3.0). There were none through 0.4.0.
• A scalar call costs more than it did. Pareto.pdf(0.5) is about 2.6us against roughly 0.2us of plain arithmetic before, because every call goes through the array dispatch. The same ten thousand points in one call take 0.3ms, about eighty times less than the loop. Scalars are carried as numpy.float64 rather than 0-d arrays and guarded with an ordinary conditional rather than np.where, which is what keeps it to 2.6us and not 5.3.
• A few seeded variates move by one unit in the last place. NumPy's vectorised log and pow round differently from its scalar ones. For Weibull(k=0.5) that is 14 values in 20,000; sample means are unchanged to twelve digits and distributions are identical. Frechet, Cauchy, Pareto, GEV_Frechet and the generalized Pareto are bit-identical.
• vectorized.accelerated answers per method, not per family. It reports True for LogNormal.pdf, where it used to report False: the density is elementary and vectorises, and only the probabilities need the error function. It also reports True for a distribution defined outside this library, because there is no longer a table of known families to consult.

Changed

The distributions. pdf, cdf, sf and ppf on all twelve continuous families accept a number, a sequence or an array and return the same kind of thing they were given. Each formula is written once, against NumPy, so the scalar and array results come from the same expression; they can still differ in the last bits, and the tests hold them to 1e-13 relative rather than to equality. rvs inverts its uniforms in a single call, which makes sampling 21x to 31x faster, and still returns a list -- returning an array would silently turn sample_a + sample_b from concatenation into elementwise addition. ppf checks its whole input before computing any of it and names an offending value.

The tail index estimators. hill_estimator is about 5x faster and hill_plot about 10x on a sample of 200,000, the plot now deriving every threshold from one cumulative sum instead of re-summing the log-excesses for each k. The two agree exactly at every k, which they did not before. trimmed_hill_plot uses suffix sums rather than subtracting from a running total. _normalised_log_spacings returns an array, and adaptive_trim_selection returns a Python float for gamma.

The multivariate families. mahalanobis, logpdf and pdf take one point or many. Over 100,000 points the quadratic form takes 5.7ms against about 2.0s a point at a time. fit_multivariate_t is about 45x faster -- 0.558s to 0.012s on 5,000 observations in three dimensions, and 4.371s to 0.103s when nu is chosen over ten candidates -- for the same answer, to ten decimals of log-likelihood. _solve_lower is gone.

The actuarial machinery. panjer_recursion evaluates its inner sum as a dot product: on a 2,000 point grid, 0.176s to 0.014s, with the probabilities unchanged to the bit on the cases measured. limited_expected_value runs its 512-node midpoint rule in one call, _severity_second_moment its 4,096-node fallback, and discretise_severity evaluates each cell edge once instead of twice. PolicyTerms.payment and LayeredSeverity.cdf, sf, ppf take one value or an array.

heavytails.vectorized is a thin shim, 447 lines to 198. It held a hand-written NumPy kernel for each of 32 (family, method) pairs, every one a transcription of the method it shadowed; the methods take arrays themselves now. The public functions behave as before, and dist.cdf(values) is the same call for new code.

Four families keep a Python loop for their probabilities. LogNormal, StudentT, InverseGamma and BetaPrime need the error function or an incomplete beta or gamma, which NumPy does not provide, so their cdf, sf and ppf cost what they always did. Their densities are elementary and are vectorised like everything else.

The package no longer describes itself as pure Python. CITATION.cff, .zenodo.json, AUTHORS.md, mkdocs.yml and the citation pages carry the title heavytails: A Python Library for Heavy-Tailed Probability Distributions; the archived 0.1.0-0.4.0 records keep the old one, which is correct, because those releases were pure Python. SECURITY.md no longer justifies the security posture with "dependency-free (pure Python)", and docs/about/license.md names NumPy's licence rather than claiming there are none to name.

Continuous integration. Coverage runs at release only. It cost about 6.4x the runtime of the suite -- 3m12s without it, 20m37s with -- and caught none of the defects found while writing this release, every one of which was in covered code. A push to main no longer re-runs the full interpreter matrix that the pull request just ran, and superseded runs are cancelled. Test jobs went from 17-22 minutes to about 3, and a merge from 19.5 minutes to 4.7.

Added

• Tests tying the citation metadata together: CITATION.cff must name the version in pyproject.toml, and the citation guidance must cite that same version. Neither was checked, which is how the documentation sat a release behind through all of 0.4.0.

Fixed

• GeneralizedPareto.sf returned exactly zero in the far tail. It computed 1 - cdf; at xi=0.5, sf(1e12) gave 0.0 where the true value is 4.0e-24. This is the tail the distribution exists to describe.
• streaming.py kept its own transcription of the Hill and moment formulas, described as matching the batch version "operation for operation" -- which held only until the batch version changed its summation order. It calls those functions now, so the agreement is structural.
• test_pickands_estimator_basic drew an unseeded sample and required a single Pickands estimate to land in (0.1, 2.0). At k=20 on 1000 points that estimator has a standard deviation of about 0.44 around a true 0.5, so the bound failed for 18% of seeds. It now averages 200 replications and asserts the estimator is centred on the true index.
• scripts/tail_index_study.py stamped its results with the installed package version, which an editable install reports from whenever its metadata was last built -- 0.2.0 against a working tree at 0.4.0. A saved study carrying the wrong version is worse than one carrying none.
• The citation guidance cited 0.3.0 through the whole 0.4.0 release, and CITATION.cff carried a placeholder for a DOI Zenodo had already minted.
• heavytails.performance evaluated densities through np.vectorize, a Python loop wearing an array interface, and heavytails.extensions guarded NumPy in the same try as SciPy, so a missing SciPy set np to None while NumPy sat installed and working.
• tests/test_zenodo_metadata.py held a third copy of the citation title rather than comparing the two files it was named for.

0.4.0 - 2026-08-22

Twelve new modules and estimators, and a dozen numerical corrections. The corrections are the reason to upgrade rather than the additions: several functions returned values that were wrong rather than merely imprecise, and one returned negative probabilities.

Some values change, because the old ones were wrong.

• GeneralizedPareto.cdf, sf and pdf return different values below mu, where they previously reported probabilities outside [0, 1].
• Weibull.pdf(0.0) returns inf for shape below one, where it previously raised ZeroDivisionError.
• InverseGamma.cdf, BetaPrime.sf, Cauchy.cdf and LogNormal.cdf return different values in the tails, where they previously returned exactly zero or kept only a handful of digits.
• cdf and sf no longer sum to exactly one for several families. Each is now computed from the branch where its own value is the small quantity, which is what lets both tails be accurate; insisting on exact complementarity would force back the subtraction that made the lower tail return zero.

Added

• heavytails.timeseries, separating heavy tails that come from volatility clustering from heavy tails that come from the shocks (#307). GARCH11 with normal or Student-t innovations, fit_garch11 by maximum likelihood, standardised residual extraction, extremal_index by the Ferro-Segers intervals estimator, and decluster by runs.

The point the module exists to make, measured on 60,000 simulated returns: Gaussian innovations -- which have no power-law tail at all -- produce a series whose estimated tail index is 4.5. Student-t innovations with 5 degrees of freedom produce 2.85. Estimating a tail index on raw returns measures the volatility process as much as the shocks, and the docstring says which of the two a given question needs.

• heavytails.copula, dependence separated from the margins (#306). GaussianCopula, StudentTCopula, GumbelCopula and GalambosCopula, with densities, exact sampling and closed-form tail dependence, plus empirical_tail_dependence for estimating it from data.

The Gaussian is included as the cautionary case: its tail dependence is zero at every correlation short of one, so fitting it to data whose extremes do arrive together reproduces the correlation faithfully and understates the joint risk with no sign that it has.

empirical_tail_dependence documents its own bias with measurements rather than a warning: on 400,000 draws from a Gaussian copula, whose true coefficient is exactly zero, it reports 0.47 at the 0.90 level and still 0.145 at 0.9995. It cannot distinguish asymptotic independence from moderate tail dependence, and says so.

Vine copulas are deliberately not included; they need pair-copula construction and structure selection and are a subsystem rather than a class.

• heavytails.multivariate, the elliptical family and joint tail dependence (#305). MultivariateStudentT and MultivariateNormal are normal scale mixtures sharing one construction, with density, sampling, Mahalanobis distance, marginals and EM fitting. tail_dependence_coefficient gives the bivariate coefficient in closed form, validated against simulation.

There is deliberately no distribution function: the multivariate t has none in closed form above one dimension, so cdf_monte_carlo estimates it and reports a standard error rather than presenting a simulation as exact. The linear algebra is pure Python, which is right at the dimensions tail dependence is asked about and would be wrong at hundreds.

• Runnable examples in the five core modules that had none (#337): 86 across heavy_tails, extra_distributions, discrete, plotting and _special. Every value was computed rather than transcribed, and tests/test_doctests.py now guards all nine modules rather than the four that previously had any.

• heavytails.vectorized, evaluating pdf, cdf, sf and ppf over many points at once using NumPy when it is installed (#308). Measured on the public call at 100,000 points, the 32 accelerated calls run 2.1x to 6.4x faster, median 3.9x. Eight families have kernels; LogNormal, StudentT, InverseGamma and BetaPrime cannot, because NumPy has neither the error function nor the incomplete beta and gamma, and accelerated() reports which is which rather than leaving a caller to guess. Without NumPy everything falls back to the loop. scripts/vectorization_benchmark.py produces the table.

• heavytails.streaming, tail index estimation over a stream without holding the sample (#310). TopK maintains the largest values in O(k) memory with a min-heap; StreamingTailIndex builds the Hill and moment estimators on it and returns bit-for-bit the batch result, since both depend on the sample only through its top k + 1 values. WindowedTailIndex does the same over the most recent observations, in O(window) memory, which is inherent rather than a gap: when the largest value in a window expires the new largest can be any of the survivors.

• tests/reference_values.json and scripts/generate_reference_values.py, a database of 209 values computed by mpmath at 50 decimal digits from the mathematical definitions, covering all twelve continuous families (#311). Every point carries its condition number, so the tolerance follows from how well the quantity can be determined at that input rather than from a flat guess. Reading the table needs no mpmath and runs in about a second.

• tests/test_distribution_properties.py, applying the generic properties -- quantile inversion, survival complementarity, monotonicity, non-negative density, support, reproducibility -- to every family from one registry rather than to four families by hand. Adds checks of the documented family relationships, of the tail index each family is supposed to have, and of moments existing exactly when the theory says.

• _gammainc_upper_reg, the regularized upper incomplete gamma computed as itself rather than as 1 - P (#309). The subtraction cannot express a result below about 1e-16, so at a=2, x=50 it returned exactly zero where the true value is 9.8e-21.

• _gammaincinv_reg, the inverse of the regularized incomplete gamma in both the lower and upper senses, so a caller can go through whichever of the two is its small quantity.

• adaptive_trimmed_hill_estimator and adaptive_trim_selection, choosing the trimming parameter for the trimmed Hill estimator from the data (#321). Each normalised log-spacing is tested against the mean of the deeper ones, with an exactly computable null distribution, and the scan trims past the deepest anomaly. The median trimming chosen equals the number of planted outliers at 0, 1, 2, 3, 5 and 8 of them, and on clean data the standard deviation is 0.0295 against 0.0292 for the plain Hill estimator. Completes the eleven estimator benchmark suite; scripts/tail_index_study.py now runs twelve.

• tests/test_doctests.py, checking that the docstring examples in the numerical modules actually reproduce. The main suite does not collect doctests, so four had stopped working unnoticed.

• heavytails.actuarial, building the aggregate loss distribution from a frequency model and a severity, and pricing the reinsurance written on it (#304). Frequency models Poisson, NegativeBinomial and Binomial -- the (a,b,0) class the Panjer recursion is defined on. PolicyTerms and LayeredSeverity apply deductibles, limits and coinsurance on either a per-loss or a per-payment basis. panjer_recursion gives the whole aggregate distribution; simulate_aggregate_loss and EmpiricalAggregate give the same interface from a sample, for the cases the recursion cannot reach. compound_moments reports the exact mean and variance, including inf when the severity is heavy enough that they do not exist. limited_expected_value, excess_of_loss_premium and AggregateLoss.stop_loss_premium price layers. Completes the actuarial item of roadmap Phase 4.

• heavytails.viz, rendering the diagnostics with matplotlib behind a new plot extra (#302): plot_tail, plot_qq, plot_hill, plot_trimmed_hill, plot_mean_residual_life and plot_parameter_stability. Each takes an optional ax and returns it, so a panel of diagnostics composes normally, and plot_tail can overlay a fitted distribution against the empirical curve. heavytails.plotting keeps returning coordinates and stays free of third-party imports, so the library itself still never requires matplotlib. Install with pip install "heavytails[plot]".

• heavytails.threshold, with mean_residual_life, parameter_stability, select_threshold and return_level (#300). Choosing the threshold dominates a peaks-over-threshold analysis and no rule settles it, so the two diagnostics come first and the automatic rule is documented as a starting point. select_threshold uses a goodness-of-fit test whose p-values are conservative, which biases it towards thresholds that are too low. return_level reports a bootstrap interval whose measured coverage is about 0.88 against a nominal 0.95, falling to 0.76 on smaller samples, because it captures sampling variability but not the error in choosing the threshold.

• heavytails.risk, with value_at_risk, expected_shortfall, tail_conditional_expectation, monte_carlo_tail_risk and mean_exists (#303). Expected shortfall has closed forms for the Pareto, log-normal, generalized Pareto and Weibull families and falls back to quadrature on the quantile function otherwise; the two paths share no code and agree to four decimal places. Expected shortfall returns inf whenever the distribution has no finite mean, rather than a large number that would look like a result, and the Monte Carlo estimator always reports standard errors.

Changed

• InverseGamma.ppf and BetaPrime.ppf invert the incomplete gamma and beta directly instead of bracketing and solving against their own distribution functions. Round-trip accuracy across the quantile range improves from 2.2e-05 to 9.7e-15 and from 5.6e-04 to 3.5e-15 respectively.

• _betaincinv_reg starts from the small-y asymptote rather than bisecting the whole exponent range, cutting StudentT.ppf from 127 to 35 microseconds with identical accuracy.

Fixed

• GeneralizedPareto returned negative probabilities below mu for every sign of xi: its validity check tested only 1 + xi z > 0, which is the upper endpoint of a bounded distribution and is satisfied far below the support. cdf(mu - 1) returned -2.586 at xi=0.4, mu=1. Found by the new vectorisation tests, and now covered by a generic property over every family.

• Weibull.pdf(0.0) raised ZeroDivisionError for shape below one, where the density diverges. It returns inf.

• Cauchy.ppf used tan(pi(u - 1/2)), whose argument sits next to +-pi/2 where the tangent is arbitrarily steep. At u = 1e-9 it returned -318309868.8 for a true -318309886.2, wrong in the eighth digit. Now uses the cotangent form, accurate to 1.9e-16.

• _phi_inverse refined its rational approximation using 0.5*(1 + erf(x/sqrt(2))), which cancels once x is a few units negative, so the correction was computed from noise. The relative error at u = 1e-12 was 4.4e-07 -- worse than the unrefined approximation. Now 4.2e-16 across the range, which also improves every normal-quantile consumer, including LogNormal.ppf and the confidence intervals in tail_index and threshold.

• GeneralizedPareto.ppf, BurrXII.ppf and Weibull.ppf formed 1 - u before taking a power or a logarithm, losing the lower tail to about seven digits. Now use log1p/expm1.

• GeneralizedPareto.cdf, BurrXII.cdf, Weibull.cdf, Cauchy.cdf and LogNormal.cdf computed a small probability as 1 - g(x) with g approaching one. Each now takes the branch where its own value is the small quantity.

• The safeguarded Newton iteration in _betaincinv_reg narrowed its bracket after computing the midpoint to fall back to, so on the first iteration -- where the starting point is the midpoint by construction -- the fallback returned that same point and the no-progress check declared convergence. I_y(50, 0.3) = 1e-3 came back with a relative error of 0.18. The fixed bisection to 1e-13 that preceded it had hidden this.

• InverseGamma.cdf computed 1 - P and returned exactly zero throughout the lower tail; cdf(0.02) at alpha=2, beta=1 is 9.8e-21, not 0.

• BetaPrime.sf computed 1 - cdf, which is exactly zero above about x=1e17 because x/(x+s) rounds to 1 there. It now uses the mirrored incomplete beta, whose argument s/(x+s) is computed rather than subtracted.

• Four docstring examples reported values their code never produced. Two claimed a tail index estimate of 0.5 on samples too small for the estimator's own sampling variability, so the figure shown had been transcribed from what the estimator should give rather than measured; both produced 0.4.

0.3.0 - 2026-08-21

Added

• bias_reduced_hill_estimator, with second_order_rho, second_order_beta and recommended_rho_k (#329). The Hill estimator's bias at large k is systematic rather than random, so it can be estimated and subtracted; with rho supplied the measured bias falls by a factor of four to eighty. second_order_rho is documented as unstable, because it is: sweeping k on a sample whose true rho is -1 gives estimates from -0.07 to -20.5, the latter at a pole. Supplying rho is strongly preferred, and the correction still helps when it is estimated.

• gpd_mle_estimator and fit_generalized_pareto (#327), the parametric peaks-over-threshold counterpart to the semiparametric estimators. Fitting uses the reduction of Grimshaw (1993), which turns the two-parameter likelihood into a one-dimensional search, so no third-party optimiser is needed; the fit matches scipy.stats.genpareto.fit to four decimal places for positive, near-zero and negative shape. Being a general-EVI estimator it handles a bounded tail, where the whole Hill family cannot.

• harmonic_moment_estimator and t_hill_estimator (#325). Hill's contributions grow without limit, so one extreme observation moves the estimate arbitrarily far; these use the bounded reciprocal ratios u / X_(i) instead. Sending a single observation of ten thousand from 1e2 to 1e30 moves the Hill estimate from 0.502 to 0.631 and moves these not at all. beta trades robustness against efficiency, and the estimator tends to the Hill estimator as beta tends to zero.

• scripts/tail_index_study.py records the heavytails version, the git commit, the Python version and the run configuration in its JSON output, so a results file can be traced back to the code that produced it. The commit is recorded as well as the version because importlib.metadata reports what is installed, which lags a working tree after a version bump.

• trimmed_hill_estimator and trimmed_hill_plot, following Bhattacharya, Kallitsis and Stoev (2019) (#323). Replacing three observations out of ten thousand with outliers moves the Hill estimate from 0.50 to 0.60; trimming five recovers 0.50. On clean data trimming ten raises the standard deviation only from 0.0296 to 0.0302.

• estimator_kwargs on tail_index_confidence_interval, so tuning parameters such as the trimming level r and the smoothing parameter u reach the estimator. Without it trimmed_hill silently ran at r = 0, which is the ordinary Hill estimator and gives no robustness at all.

• Contaminated scenarios in scripts/tail_index_study.py, since robustness is invisible on clean data.

Changed

• The tail_index module docstring now states that every estimator returns the extreme-value index gamma = 1/alpha, not the tail index alpha (#322). The module name refers to the quantity estimated, not the parameterisation returned, and the two conventions are reciprocals.

• smoothed_hill_estimator, the smoothed Hill estimator of Resnick and Stărică (1997), and smoothed_hill_variance_ratio, which reports the asymptotic variance reduction it achieves (#319). Averaging the Hill estimate over j in (k, u*k] reduces its asymptotic variance by 39% at u = 2 and 55% at u = 3, and in the simulation study it has the lowest RMSE of any estimator here on a heavy tail. It inherits Hill's restriction to positive gamma.

• generalized_hill_estimator, the UH estimator of Beirlant, Vynckier and Teugels (1996) (#299). Unlike the Hill estimator it is consistent for every extreme-value index, not only positive ones. On a Uniform(0,1) sample, whose index is -1, Hill can only ever return a positive number and reports about +0.026; the generalized Hill estimator recovers -0.99.

• hill_plot, which sweeps k on a logarithmic grid and returns the (k, gamma) series. The documentation already told readers to find a plateau; they now have something that produces one.

• tail_index_confidence_interval, with an asymptotic interval for the Hill estimator and a percentile bootstrap for all four. Requesting the asymptotic interval for an estimator that has no established closed form raises rather than reporting a number with no basis.

• scripts/tail_index_study.py, a simulation study reporting bias, standard deviation and RMSE for every estimator across known indices and sample sizes. Its results are summarised in the tail estimation guide.

• _phi_inverse moved from heavy_tails to _special, alongside the other numeric helpers, so tail_index can use it without importing the distribution module. It was never public.

0.2.0 - 2026-08-21

Added

• scripts/special_function_accuracy.py, which sweeps both special functions against mpmath at 50 decimal digits and reports the worst relative error. Both are accurate to about 12 significant digits over the ranges the distributions use, with no degradation at the series/continued-fraction switch. tests/test_special_accuracy.py asserts those bounds.

• Kolmogorov-Smirnov and Anderson-Darling goodness-of-fit tests in heavytails.validation.GoodnessOfFitTests, which previously raised NotImplementedError (#301). Both are reported by AutoFit.compare_distributions and by the heavytails compare command, so a comparison now says whether the winning family fits, not only how it ranks.

• ConvergenceError, raised when a solver cannot reach its tolerance. It previously returned its best guess, which the caller could not distinguish from a converged answer.

• LogNormal.ppf raised OverflowError when the quantile exceeded the float range instead of returning inf. The median of LogNormal(mu=1000) is exp(1000), which is genuinely not representable, so inf is the correct answer; raising broke parameter sweeps and made the failure look like a caller error. (#296)

• LogNormal.sf is computed with math.erfc rather than as 1 - cdf(x), which reached exactly zero by x = 1e5 and carried no information beyond it. It now matches SciPy to about 1e-14 relative out to x = 1e12.

• Python 3.13 support, covered by the CI test matrix and declared in the package classifiers.

• py.typed marker, so the type annotations that ship with the package are visible to downstream type checkers. The Typing :: Typed classifier was previously advertised without one.

• heavytails.__version__, resolved from the installed distribution metadata.

• --version / -V flag on the heavytails command-line interface.

• cli installation extra (pip install "heavytails[cli]"). The console script depends on typer and rich, which were previously development-only dependencies, so the entry point was broken for anyone installing from PyPI.

• CodeQL analysis and dependency-review workflows.

• .github/CODEOWNERS, .gitattributes and .zenodo.json.

• StudentT.cdf, StudentT.sf and StudentT.ppf. The class previously offered only pdf and rvs, with a docstring stating that the CDF and PPF "require special functions not in stdlib" — but the regularized incomplete beta needed to write them was already implemented in extra_distributions. All three agree with SciPy to around 1e-14.

• Cauchy.sf, Frechet.sf and GEV_Frechet.sf, so every continuous family now provides the full interface the documentation advertises.

• YuleSimon.sf and YuleSimon.ppf. The survival function uses the closed form P(X > k) = k * B(k, rho + 1), and the quantile function brackets and bisects rather than scanning linearly.

• heavytails._special, holding the shared numeric special functions so that both distribution modules can use them without an import cycle. The previous names remain importable from heavytails.extra_distributions.

• _betaincinv_reg, an inverse for the regularized incomplete beta. It solves in log-space with a symmetry reduction, which is what keeps extreme quantiles accurate.

• Documentation pages that existing pages already linked to but which had never been written: CLI reference, diagnostics guide, extreme value theory, tail index estimation theory, validation studies, architecture, benchmarking, code review, and an executable basic_usage notebook.

Changed

• Migrated project metadata to the PEP 621 [project] table.

• Grouped Dependabot updates so routine bumps arrive as a few reviewable pull requests rather than one per package.

• Continuous integration now also runs on pull requests targeting develop, builds and metadata-checks the distributions before publishing, verifies that poetry.lock matches pyproject.toml, and builds the documentation with --strict.

• Security scanning fails the build on findings instead of uploading a report that no one reads. safety was replaced by pip-audit, which needs no account to run.

• Refreshed the locked dependency set, clearing 116 known vulnerabilities reported against the previously locked development and documentation toolchain.

• Pre-commit hooks are pinned to the same tool versions as the development dependencies, so local hooks and CI now agree.

• Replaced the Makefile targets, which measured coverage of scripts/ rather than of the package.

• mkdocstrings is configured for Google-style docstrings, which is what the package actually uses. Under the previous numpy setting no Args: or Returns: section was parsed anywhere in the API reference.

• The documentation navigation lists every page. Five API reference entries all pointed at the same directory, and the gen-files script duplicated the hand-written reference pages while literate-nav looked for a SUMMARY.md that was never generated.

Fixed

• The continued-fraction branch of the regularized lower incomplete gamma was wrong, not merely imprecise: its Lentz recurrence was missing its leading term and its b was shifted by one, so P(20, 21) returned 0.0 against a true 0.6157 (#297). That branch is reached whenever x >= a + 1, and InverseGamma.cdf evaluates P(alpha, beta/x), so its reported probabilities were wrong by factors of 2 to 17 in the lower tail. Found by comparing against mpmath; the property-based checks could not have caught it, because a consistently wrong value is still monotone and still lies in [0, 1].

• Quantile functions now behave the same way across every family (#298). Previously the answer at the edges depended on which family you asked: Pareto, Weibull, Frechet, GEV_Frechet, BurrXII and LogLogistic raised OverflowError where the quantile simply exceeded the float range, InverseGamma raised a ValueError about failing to bracket the root, LogNormal returned inf, and BetaPrime returned a finite value. All of them now return inf.

• DiscretePareto.ppf accepted any float, including 0.0, 1.5 and negatives, and silently returned the support bounds. It now raises ValueError like every other family.

• _ppf_monotone narrowed its bracket only on bisection fallbacks, never on accepted Newton steps, so a run of Newton steps could exhaust the iteration budget with the bracket as wide as it started. It now narrows on every iteration, which also fixed a non-monotonicity in InverseGamma.ppf.

• heavytails benchmark raised ZeroDivisionError on platforms with a low-resolution wall clock, because a sub-millisecond timing measured exactly zero seconds. Timings now use time.perf_counter().

• The same defect in the performance tests, which failed intermittently on Windows.

• Removed an invalid PyPI classifier (Topic :: Scientific/Engineering :: Statistics) that would have been rejected on upload.

• scripts/pyproject_updater.py depends on tomlkit and packaging, which were never declared and only happened to be installed transitively.

• Pinned Poetry 2.2.1 in every workflow. The previous 1.8.3 pin cannot read the version 2.1 lock file this repository uses.

• Removed a preferred-citation entry from CITATION.cff that pointed at an unpublished paper with a placeholder DOI, which citation tooling would have emitted as a real reference.

• YuleSimon.pmf raised OverflowError for k of about 170 and above, because it multiplied gamma functions that overflow individually even though their ratio is small. Since sampling called it in a loop, drawing from the tail crashed. It is now evaluated with lgamma.

• Nineteen documentation links pointed at pages that did not exist. The documentation job now builds with --strict, so a broken link fails CI.

• convergence_validation took a _max_iter parameter that was never used and documented it under a different name.

• Removed a stale xfail marker on the Student-t PPF convergence test. It was recording precision loss that the new incomplete beta inverse eliminates.

Removed

• roadmap.improved_incomplete_beta, a placeholder that delegated to _betainc_reg unchanged. The accuracy work it stood in for is now done and measured.

• roadmap.safe_lognormal_ppf, a workaround that caught the overflow above. It is redundant now that LogNormal.ppf handles the case itself.

• tox.ini, which ran the test suite against scripts/ instead of the package and duplicated the CI matrix.

• IMPROVEMENT_PLAN.md from the repository root. It described gaps that have since been closed; forward-looking plans live in ROADMAP.md.

• The isort development dependency and its configuration, superseded by Ruff's I rules.

0.1.0 - 2025-10-25

Added

• Continuous heavy-tailed distributions implemented from first principles: Pareto, Cauchy, Student-t, Log-Normal, Weibull, Fréchet and GEV (ξ > 0).
• Additional continuous families: Generalized Pareto, Burr XII, Log-Logistic (Fisk), Inverse-Gamma and Beta-Prime.
• Discrete heavy-tailed distributions: Zipf, Yule–Simon and Discrete Pareto.
• Tail index estimators: Hill, Pickands and moment.
• Diagnostic plotting helpers for log–log tail and QQ plots.
• Deterministic RNG wrapper for reproducible sampling.
• Custom incomplete-gamma and incomplete-beta implementations, and a safeguarded-Newton numeric PPF solver for families without a closed form.
• heavytails command-line interface.
• Documentation site built with MkDocs Material.