diff --git a/.github/workflows/actions.yml b/.github/workflows/actions.yml index 65be7595..159bce14 100644 --- a/.github/workflows/actions.yml +++ b/.github/workflows/actions.yml @@ -6,6 +6,11 @@ run-name: SurPyval CI # master (i.e. release merges + tags). Feature work accumulates on `develop` # via PRs; day-to-day pushes to feature branches no longer each trigger the # full suite. See docs/Contributing.rst for the branching / release flow. +# +# The jobs below are not all gated the same way. Lint runs on everything; +# the test suite and the documentation build run only on the release pull +# request into `master`, where the time they cost is worth paying. Each +# job carries the reasoning for its own condition. on: pull_request: branches: [develop, master] @@ -40,7 +45,26 @@ jobs: - name: black run: black --check $SRC + # The test suite runs on the release pull request into `master`, and on + # pushes to `master`, but not on pull requests into `develop`. + # + # It is roughly nine minutes across the three interpreters, against + # about one for lint, and paying that on every feature pull request + # made the edit-review loop the slowest part of working on the + # package. Development is single-maintainer and the suite is run + # locally before pushing, so the pull-request run was mostly + # confirming what was already known. + # + # What that gives up is real, and worth naming: a failure that only + # appears on one interpreter is now found when the release is + # prepared, with a release's worth of commits to search rather than + # one. That is not hypothetical -- the doctest comparison landed + # green on 3.11 and broke on 3.12 and 3.13, and it was the + # pull-request run that caught it. Run the suite locally across more + # than one interpreter when touching numerics, or open the release + # pull request early and let it sit. surpyval_ci: + if: github.event_name != 'pull_request' || github.base_ref == 'master' runs-on: ubuntu-latest strategy: fail-fast: false @@ -82,6 +106,36 @@ jobs: python -m pytest -n auto --cov=surpyval --cov-report= --ignore=surpyval/tests/alpha --run-ml + # Execute the ``>>>`` examples in the docstrings and compare their + # printed output. This is what a user (or an agent) sees from + # ``help(Weibull.fit)``, so it is a documented promise like any + # other, and it went stale silently until it was checked. Kept as + # its own step -- separate from the suite above -- so a failure + # reads as "the docs drifted", not "a test broke", and so the + # doctest collection does not disturb coverage or xdist. + # + # ``surpyval/tests`` is excluded because the test modules have no + # user-facing examples, and ``surpyval/alpha`` because it is not + # part of the release contract (same reason the suite skips it). + # Option flags are set in pyproject.toml; conftest.py adds the + # numeric comparison that lets the examples record real output. + - name: doctests + run: + python -m pytest --doctest-modules surpyval + --ignore=surpyval/tests --ignore=surpyval/alpha + + # The numeric comparison above only runs for an example whose + # output has actually drifted -- a handful on any one machine, and + # a different handful on each. That leaves a gap in it invisible + # until CI hits the one example that needed it. This second pass + # routes every example through it, so the fallback is checked + # against all of them. Fifteen seconds. + - name: doctests (numeric comparison forced) + run: + python -m pytest --doctest-modules surpyval + --ignore=surpyval/tests --ignore=surpyval/alpha + --doctest-force-numeric + - name: coverage run: | coverage report @@ -93,3 +147,54 @@ jobs: with: name: coverage-html-report path: htmlcov/ + + # Build the documentation on the release pull request only. + # + # The docs execute every ``.. jupyter-execute::`` cell during the + # build, so they are a second test suite that runs the public API for + # real -- and one that can be broken by a change that touches no + # documentation file at all. Removing Gamma's probability-plot fitting + # did exactly that: a cell looping over the fit methods started + # raising, and nothing noticed, because Read the Docs builds only + # ``master`` and tags. The failure would have appeared as a broken + # hosted build after the release rather than on the change that caused + # it. + # + # ``github.base_ref`` is set only for pull_request events, so this + # runs on the develop -> master release pull request and nowhere else. + # It is deliberately not run on pushes to master: Read the Docs + # rebuilds there anyway, and by that point the gate has nothing left + # to gate. + # + # The point is to reproduce the hosted build, so this matches + # .readthedocs.yaml rather than the test jobs above: Python 3.12, and + # the package installed with its own ``docs`` extra. + docs: + if: github.event_name == 'pull_request' && github.base_ref == 'master' + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v4 + + - name: Set-up Python + uses: actions/setup-python@v5 + with: + python-version: '3.12' + + - name: Install dependencies + run: | + python -m pip install --upgrade pip + pip install ".[docs]" + + # Not -W: there are pre-existing warnings (duplicate changelog + # labels from autosectionlabel, an rtd-theme deprecation, three + # autodoc imports). Clearing those is worth doing, and then + # failing on warning here and in .readthedocs.yaml together -- + # turning it on before then would just fail every release. + - name: sphinx build + run: python -m sphinx -b html docs docs/_build/html + + - name: Upload documentation artifact + uses: actions/upload-artifact@v4 + with: + name: docs-html + path: docs/_build/html/ diff --git a/.gitignore b/.gitignore index 9a75e36c..f28a9153 100755 --- a/.gitignore +++ b/.gitignore @@ -24,3 +24,6 @@ __pycache__/ publish.sh + +# Local multi-interpreter check environments (scripts/check_all_pythons.py) +.venvs/ diff --git a/.readthedocs.yaml b/.readthedocs.yaml index 3b8c8153..c014b91d 100644 --- a/.readthedocs.yaml +++ b/.readthedocs.yaml @@ -34,4 +34,5 @@ python: install: - method: pip path: . - - requirements: docs/requirements.txt \ No newline at end of file + extra_requirements: + - docs \ No newline at end of file diff --git a/conftest.py b/conftest.py index 1df740b0..227d5e0b 100644 --- a/conftest.py +++ b/conftest.py @@ -1,4 +1,8 @@ -"""Opt-in gating for the slow parts of the test suite. +"""Suite-wide fixtures: doctest number comparison, and opt-in gating. + +The first half of this file makes the ``--doctest-modules`` run compare +the numbers in an example's output as numbers rather than as text; see +the comment above ``RTOL``. The rest is the opt-in gating below. Two groups are skipped unless asked for, because both are expensive and neither guards a regression that the default run would miss quickly: @@ -31,8 +35,127 @@ the run dies on "unrecognized arguments". """ +import doctest +import math +import re + import pytest +# --------------------------------------------------------------------------- +# Numeric comparison for the ``--doctest-modules`` run +# --------------------------------------------------------------------------- +# doctest compares printed output as text. That is the wrong test for a +# library whose examples end in an optimiser: the same fit lands on +# ``b = 4.1995e-05`` under one Python and ``4.2032e-05`` under the next, +# and numpy prints eight significant digits either way, so a byte-exact +# comparison fails on a difference no reader would call a difference. +# +# The alternative -- trimming every documented number to the digits that +# happen to agree everywhere -- makes the docstring show something the +# user's own session will not produce, which is the thing these examples +# exist to avoid. So the examples record the real output, in full, and +# the numbers in it are compared as numbers. +# +# The fallback only runs after the ordinary text comparison has failed, +# and only fires when the two outputs are identical apart from their +# numeric literals -- same words, same brackets, same integer-vs-float +# shape ("1" never matches "1.", which is a dtype change worth +# failing on). What it forgives is the value drifting inside a +# tolerance. What it still catches is everything that actually went +# wrong when this was first switched on: a stale value from another +# parameterisation, a different function being called, the wrong array +# shape, an exception, a missing import. +# +# RTOL is set by the loosest genuine disagreement between supported +# Pythons -- the Duane example above, at 9e-4 -- with no margin beyond +# that. ATOL exists for the one other case, a restoration factor whose +# true value is zero and which surfaces as 1e-16 with whatever sign and +# mantissa the optimiser stopped on; relative tolerance is meaningless +# there. +RTOL = 1e-3 +ATOL = 1e-12 + +_NUMBER = re.compile(r"[-+]?(?:\d+\.\d*|\.\d+|\d+)(?:[eE][-+]?\d+)?") +_BLANKLINE = re.compile(r"(?m)^%s\s*?$" % re.escape(doctest.BLANKLINE_MARKER)) + + +def _skeleton(text: str) -> str: + """The text with each number replaced by its *kind*. + + Integers and floats get different placeholders so that a change in + dtype -- ``array([1, 2])`` becoming ``array([1., 2.])`` -- is still + a failure rather than two numbers that happen to be equal. + + Whitespace is dropped entirely. numpy pads an array's columns to its + widest element, so shortening one number moves the spaces around + every other: ``[ 6.32508961 17.37701969]`` against + ``[ 6.3250866 17.377018 ]``. Those spaces carry no meaning the + numeric comparison below has not already made. + """ + + def mark(match: re.Match) -> str: + token = match.group(0) + return "~f" if ("." in token or "e" in token or "E" in token) else "~i" + + return "".join(_NUMBER.sub(mark, text).split()) + + +def _numerically_equal(want: str, got: str) -> bool: + # ```` stands for an empty line in the expected output. + # The text comparison substitutes it before matching, so this one has + # to as well, or a model repr with a blank line in it can never reach + # the numeric comparison at all. + want = _BLANKLINE.sub("", want) + + if _skeleton(want) != _skeleton(got): + return False + wants = _NUMBER.findall(want) + gots = _NUMBER.findall(got) + if not wants or len(wants) != len(gots): + return False + return all( + math.isclose(float(w), float(g), rel_tol=RTOL, abs_tol=ATOL) + for w, g in zip(wants, gots) + ) + + +_text_check_output = doctest.OutputChecker.check_output + + +def _check_output(self, want, got, optionflags): + if _text_check_output(self, want, got, optionflags): + return True + return _numerically_equal(want, got) + + +# Patched on the base class rather than installed as a checker: pytest +# builds its own ``LiteralsOutputChecker`` subclass and calls up to this +# method, so overriding here survives both plain ``doctest`` and pytest, +# and does not depend on pytest's internals. +_patched = _check_output # type: ignore[assignment] +doctest.OutputChecker.check_output = _patched # type: ignore[method-assign] + + +def _forced_check_output(self, want, got, optionflags): + """As above, but the numeric path is the *only* path. + + The fallback normally runs only when an example's output has + actually drifted, which on any one machine is a handful of them. A + gap in it -- the ```` markers it did not strip, say -- + therefore stays invisible locally and surfaces in CI, on whichever + Python happens to compute a different last digit. + + Under ``--doctest-force-numeric`` every example whose output + contains a number is compared numerically instead, so the fallback + is exercised against all 229 of them rather than against today's + accidental few. Outputs with no numbers keep the text comparison; + there is nothing in them for this to compare. + """ + if not _NUMBER.search(want): + return _text_check_output(self, want, got, optionflags) + return _numerically_equal(want, got) + + OPT_IN = { "ml": ( "--run-ml", @@ -55,6 +178,16 @@ def pytest_addoption(parser): default=False, help=f"run the {description} (skipped by default)", ) + parser.addoption( + "--doctest-force-numeric", + action="store_true", + default=False, + help=( + "compare every doctest example's numbers numerically, not " + "only those whose text has drifted; exercises the fallback " + "against all of them" + ), + ) def pytest_configure(config): @@ -62,6 +195,10 @@ def pytest_configure(config): config.addinivalue_line( "markers", f"{mark}: {description}; opt in with {flag}" ) + if config.getoption("--doctest-force-numeric"): + doctest.OutputChecker.check_output = ( # type: ignore[method-assign] + _forced_check_output # type: ignore[assignment] + ) def pytest_collection_modifyitems(config, items): diff --git a/docs/Contributing.rst b/docs/Contributing.rst index e39cc43a..3fede8b8 100644 --- a/docs/Contributing.rst +++ b/docs/Contributing.rst @@ -22,10 +22,55 @@ documentation build from running on every change: Continuous integration (``.github/workflows/actions.yml``) therefore runs on **pull requests into develop or master** and on **pushes to master**, rather -than on every push to every branch. Read the Docs is configured to build -``master`` and tags only. The net effect is that the full test suite and the -documentation build run once per pull request and once per release, instead of -once per intermediate commit. +than on every push to every branch. Not every job runs on every event: + +.. list-table:: + :header-rows: 1 + :widths: 30 70 + + * - Event + - Jobs + * - Pull request into ``develop`` + - lint only (about a minute) + * - Pull request into ``master`` (the release) + - lint, the test suite across three interpreters, and the + documentation build (about ten minutes) + * - Push to ``master`` / tag + - lint and the test suite; Read the Docs rebuilds the hosted + documentation + +The test suite and the documentation build are both gated at the release +rather than on every pull request because of what they cost: the suite is +roughly nine minutes across the three interpreters and the documentation build +around three from cold, against about one for lint. Paying that on every +feature pull request made the edit-review loop the slowest part of working on +the package, and with a single maintainer running the suite locally before +pushing, the pull-request run was mostly confirming what was already known. + +The trade-off is real and worth understanding before you rely on it. A failure +that appears on only one interpreter, or a change that breaks a documentation +example, is now found when the release pull request is opened -- with a +release's worth of commits to search through rather than one. So: + +* Run the suite locally before pushing, and across more than one interpreter + when you have touched anything numerical. ``scripts/check_all_pythons.py`` + does exactly that -- it runs what continuous integration would have run, on + 3.11, 3.12 and 3.13: + + .. code-block:: bash + + python scripts/check_all_pythons.py # all three + python scripts/check_all_pythons.py 3.12 # just one + python scripts/check_all_pythons.py --skip-install # reuse as-is + + It keeps its environments in ``.venvs/`` (git-ignored) and reuses them, so + only the first run pays for the installs. It uses ``uv`` when that is + available and falls back to ``venv`` and ``pip`` when it is not. + +* Build the documentation locally when you change the behaviour of a public + function, since documentation cells call the real API. +* On a long-running branch, open the release pull request early and let it sit, + so the full run has somewhere to fail before the release itself. Documentation ------------- @@ -41,8 +86,7 @@ To build the documentation locally: .. code-block:: bash - pip install -e . - pip install -r docs/requirements.txt + pip install -e ".[docs]" sphinx-build -b html docs docs/_build/html When writing documentation, prefer ``.. jupyter-execute::`` over static diff --git a/docs/Parametric SurPyval Modelling.rst b/docs/Parametric SurPyval Modelling.rst index 895632c7..26b3f134 100644 --- a/docs/Parametric SurPyval Modelling.rst +++ b/docs/Parametric SurPyval Modelling.rst @@ -393,19 +393,29 @@ In practice that sensitivity is the estimator's problem to solve, not yours. Shi x = surv.Gamma.random(10_000, 3.0, 2.0) + 10.0 print('Truth : gamma=10.000, alpha=3.000, beta=2.000') - for how in ['MPP', 'MOM', 'MSE', 'MPS', 'MLE']: + for how in ['MOM', 'MSE', 'MPS', 'MLE']: m = surv.Gamma.fit(x, offset=True, how=how) print('{:6s}: gamma={:.3f}, alpha={:.3f}, beta={:.3f}'.format( how, m.gamma, *m.params)) -This was not always so. Earlier versions could land on an absurd tuple - a negative ``gamma`` with a shape parameter inflated by two orders of magnitude - which was nonetheless an acceptable *distribution*, precisely because of the flat trade-off described above. Two separate causes were at work, and both are now fixed. The probability-plotting search was stranded at a poor local optimum by a single starting shape, so it now tries several and keeps the best correlation. The moment-based initialisers took their moments from the *unshifted* data, where the offset dominates every moment and the shape estimate explodes - 649 for a true shape of 3 - so they now estimate the threshold first and read the remaining parameters off ``x - gamma``. +``MPP`` is absent from that list because the Gamma does not offer it. A +probability plot needs a straight-line y-axis that can be drawn *before* +the parameters are known; the Gamma's CDF is the regularised incomplete +gamma function, with the shape inside the special function rather than +outside as an exponent, so the only such axis is the inverse incomplete +gamma — which needs the shape. To draw the axis you need the answer. +``Gamma.fit(x, how="MPP")`` raises rather than guessing a shape to draw +the axis with; ``Gamma.plot()`` is unaffected, since by then the fitted +parameters are in hand. + +This was not always so. Earlier versions could land on an absurd tuple - a negative ``gamma`` with a shape parameter inflated by two orders of magnitude - which was nonetheless an acceptable *distribution*, precisely because of the flat trade-off described above. The moment-based initialisers took their moments from the *unshifted* data, where the offset dominates every moment and the shape estimate explodes - 649 for a true shape of 3 - so they now estimate the threshold first and read the remaining parameters off ``x - gamma``. The underlying caution still stands, though, and it is worth keeping in mind for your own data: - **Judge an offset fit by what it predicts, not only by the printed parameters.** Plot it against the non-parametric estimate, or compare the survival function, quantiles, mean and variance. Two parameter tuples that look very different can imply nearly the same distribution. - **If you need ``gamma`` itself to be meaningful** - you are interpreting it as a guaranteed minimum life, say - prefer ``MLE``, which remains the most accurate on the parameters, and treat a single point estimate of a threshold with care regardless of method. -``test_offset_divergence.py`` in the test suite pins this down with measured KL and Wasserstein distances alongside parameter tolerances: ``MLE`` is held to 5% on every parameter, ``MOM`` to 10% and ``MPP`` to 20%, with the implied distributions essentially identical in all three cases. +``test_offset_divergence.py`` in the test suite pins this down with measured KL and Wasserstein distances alongside parameter tolerances: ``MLE`` is held to 5% on every parameter across the offsettable distributions, and ``MOM`` to 10%, with the implied distributions essentially identical either way. Fixing parameters ----------------- diff --git a/docs/changelog.rst b/docs/changelog.rst index 3b009543..468aeb9e 100644 --- a/docs/changelog.rst +++ b/docs/changelog.rst @@ -1,6 +1,309 @@ Changelog ========= +v0.19.1 (unreleased) +-------------------- + +- **``scripts/check_all_pythons.py`` runs the CI checks locally on every + supported interpreter.** With the suite no longer running on pull + requests into ``develop`` (below), this is the other half of the + trade: one command runs the test suite and both doctest passes on + 3.11, 3.12 and 3.13, and refuses to say "passed" unless all of them + did. + + It keeps its environments in a git-ignored ``.venvs/`` and reuses + them, so only the first run pays for the installs; it uses ``uv`` + when available and falls back to ``venv`` and ``pip`` when not, and + reports an interpreter that is not installed rather than failing on + it. The command list is deliberately a copy of the workflow's, so + what it runs is what CI would have run. + +- **The test suite runs on the release pull request, not on every one.** + Pull requests into ``develop`` now run lint only, about a minute + against the nine the suite takes across three interpreters. The suite + still runs in full on the release pull request into ``master`` and on + pushes to ``master``. + + The reason is the edit-review loop: with a single maintainer running + the suite locally before pushing, the pull-request run was mostly + confirming what was already known, and it was the slowest part of + working on the package. + + What this gives up is stated rather than glossed: a failure that + appears on only one interpreter is now found when the release is + prepared, with a release's worth of commits to search rather than one. + That is not hypothetical -- the doctest numeric comparison two entries + below landed green on 3.11 and failed on 3.12 and 3.13, and it was the + pull-request run that caught it. ``Contributing.rst`` now says which + jobs run on which event, and what to run locally to compensate. + +- **The documentation build runs in CI on the release pull request.** + The docs execute every ``.. jupyter-execute::`` cell as they build, so + they are a second test suite that exercises the public API for real -- + and one that a change touching no documentation file at all can break, + as the Gamma entry below did. Read the Docs builds only ``master`` and + tags, so until now that break would have surfaced as a failed hosted + build *after* a release. + + The new job is conditioned on ``github.base_ref == 'master'``, which + is set only for pull requests, so it runs on the ``develop`` -> + ``master`` release pull request and nowhere else. It is not run on + pushes to ``master`` either: Read the Docs rebuilds there anyway, and + by then the gate has nothing left to gate. It matches + ``.readthedocs.yaml`` rather than the test jobs -- Python 3.12, the + package installed via its own ``docs`` extra -- because its purpose is + to reproduce the hosted build, and it uploads the rendered HTML as an + artifact. + + It does not build with ``-W``; the build currently emits 18 warnings, + mostly duplicate labels from ``autosectionlabel`` meeting the + changelog's repeated section headings. Clearing those and then failing + on warning here and in ``.readthedocs.yaml`` together is worth doing + separately. + + The residual gap is deliberate: a documentation break introduced on a + pull request into ``develop`` is caught when the release is prepared, + not when it lands. Building on every pull request would cost minutes + on each, and a path filter would not have helped here -- the change + that broke the build was in ``gamma.py``, not under ``docs/``. + +- **Fixed a documentation build broken by the Gamma MPP removal.** The + offset-threshold section of *Parametric SurPyval Modelling* ran a + ``jupyter-execute`` cell looping over + ``['MPP', 'MOM', 'MSE', 'MPS', 'MLE']`` for a shifted Gamma. Since + ``Gamma.fit(how="MPP")`` now raises, that cell raised, and because + documentation cells are executed during the build the whole build + failed. + + Nothing caught it: continuous integration does not build the + documentation, and Read the Docs builds only ``master`` and tags, so + it would have surfaced as a failed hosted build at the next release + rather than on the change that caused it. It was found by running a + build to validate the ``docs`` extra below. + + The prose around the cell had gone stale the same way -- it described + the multi-start probability-plotting search that the removal deleted, + and quoted an ``MPP`` tolerance from ``test_offset_divergence.py`` + that no longer exists. It now explains why the Gamma has no + probability plot at all: the shape sits inside the regularised + incomplete gamma rather than outside as an exponent, so the only + straight-line axis is the inverse incomplete gamma, which needs the + very shape being estimated. ``Gamma.plot()`` is unaffected, since by + then the parameters are known. + +- **The documentation toolchain is a ``docs`` extra.** + ``pip install -e ".[docs]"`` now installs everything needed to build + the documentation, alongside the ``tests`` extra that was already + there. ``docs/requirements.txt`` is gone: its pins moved into + ``pyproject.toml`` verbatim, and Read the Docs installs the extra + directly via ``extra_requirements``. Keeping both would have meant two + copies of the same pinned toolchain, which is the arrangement that + drifts. + + The pins are unchanged, including the ``ipykernel==6.31.0`` cap and + the reason for it -- jupyter-sphinx notebook execution dies against + the ipykernel 7 line. ``matplotlib`` is not repeated in the extra; it + is a runtime dependency of the package, which is installed alongside. + + Part of #141. + +- **CI now runs the docstring examples.** ``pytest --doctest-modules`` + over the package is a new step in the deployment workflow, and every + one of the 229 docstring examples passes. It was 59 failing tests + when the flag was first turned on. + + A docstring example is a promise about what the library prints, and + it is the one users and coding agents reach for first -- + ``help(Weibull.fit)`` is faster than opening the docs. Nothing was + checking it, so it drifted: examples recorded the output of an + optimiser two rewrites ago, of numpy 1's scalar repr, of a module + that has since moved. + + What the run found, beyond the cosmetic drift: + + - Twelve examples could not run at all. Six regression docstrings + (``PH``, ``AH``, ``PO``, ``AFT``, ``AcceleratedLife``, ``Frailty``) + were sketches -- ``model = PH(Weibull).fit(x, Z=covariates, c=c)`` + with ``x``, ``covariates`` and ``c`` never defined. Four more used + ``>>>`` on the continuation lines of a multi-line call, so pasting + them raised ``SyntaxError``. ``plotting_positions`` imported from + ``surpyval.nonparametric``, which moved to + ``surpyval.univariate.nonparametric`` several releases ago. + ``ParametricFitter.fit`` demonstrated ``how='MPP'`` on + interval-censored input, which now (correctly) requires the Turnbull + heuristic and raises without it. All are now runnable, with data. + + - The five ``ParametricRegressionModel`` prediction examples + (``sf``, ``ff``, ``df``, ``hf``, ``Hf``) had been copied from the + univariate ``Parametric`` class and never adapted: they built a + ``Weibull.from_params([10, 3])`` and called it with no covariates at + all, documenting a signature the method does not have. They now fit + a ``WeibullPH`` and pass ``Z``. + + - ``Parametric.var()`` claimed 11.229 for a Weibull(10, 3). The + variance is 10.533 (``100 Gamma(5/3) - (10 Gamma(4/3))^2``); the + code was right. + + - Several examples fitted unseeded random data and then recorded + specific digits, which cannot be reproducible. They now seed. + + ``Parametric.hf`` and ``Parametric.Hf`` returned a 0-d array + (``array(0.012)``) for scalar input where ``sf``, ``ff``, ``df`` and + ``qf`` all returned a numpy scalar, and ``cs`` did the same; their + own ``Returns`` sections promised "the scalar value ... if a scalar + was passed". That is now true. The 0-d array came from ``np.where``, + which does not collapse. + + **The numbers in the examples are compared as numbers.** doctest + compares printed output as text, which is the wrong test for a library + whose examples end in an optimiser: the same ``Duane`` fit lands on + ``b = 4.1995e-05`` under Python 3.11 and ``4.2032e-05`` under 3.12, + and numpy prints eight significant digits either way. Sixteen of the + 229 examples disagree between those two Pythons somewhere in their + digits. + + The obvious workaround -- trimming each documented number back to the + digits that agree everywhere -- makes the docstring show something the + reader's own session will not produce, which is precisely what these + examples exist to avoid. So the examples record the real output, in + full, and ``conftest.py`` installs a fallback comparison that runs + only after the ordinary text comparison has failed. It fires when the + two outputs are identical apart from their numeric literals -- same + words, same brackets, same integer-versus-float shape, so ``1`` never + matches ``1.`` and a dtype change is still a failure -- and then + compares the numbers with ``rel_tol=1e-3``, set by the loosest genuine + disagreement between supported Pythons with no margin beyond it, and + ``abs_tol=1e-12`` for a restoration factor whose true value is zero + and which surfaces as ``1e-16`` with whatever mantissa the optimiser + stopped on. + + What that forgives is a value drifting inside the tolerance. What it + still catches is every defect listed above: a stale value from another + parameterisation, the wrong function being called, the wrong shape, an + exception, a missing import. ``surpyval/tests/test_doctest_checker.py`` + pins both halves of that, using the real output pairs observed on + different Pythons. + + The fallback only runs when an example has actually drifted, which on + any one machine is a handful of them -- and a different handful on + each. A gap in it is therefore invisible locally and surfaces in CI, + on whichever Python computed a different last digit. So the doctest + step runs twice: once normally, and once under + ``--doctest-force-numeric``, which routes every example whose output + contains a number through the numeric comparison. Fifteen seconds, and + the fallback is exercised against all 229 examples rather than + today's accidental few. + + ``NORMALIZE_WHITESPACE`` is set in ``pyproject.toml`` for the same + reason: numpy picks its own line breaks and column padding for an + array and both move with the width of the widest element. + + This closes #158. + +- **The distribution docstring examples now show what you actually see.** + ``pytest --doctest-modules`` over ``distributions/`` is green: 139 + examples, no failures. Previously 39 failed. + + Most were the numpy 2 scalar repr. ``Weibull.mean(3, 4)`` prints + ``np.float64(2.7192074311664314)``, where the docstring recorded the + bare ``2.7192074311664314`` that numpy 1 used to print. The examples + now record the wrapper, because that is what appears at a prompt -- + the alternative was ``np.set_printoptions(legacy="1.25")`` in a test + fixture, which would have kept the docstrings prettier by showing + readers something their own session will not produce. + + Four ``qf`` examples printed wider than the 79-column limit once the + real output was recorded, numpy having rewrapped the arrays. Rather + than hand-wrap them into something numpy would not emit, those + examples take fewer probabilities: what is shown is exactly what that + input produces. + + Two scalar examples had drifted in the last digit, and are re-recorded. + + With the examples now true, ``--doctest-modules`` is worth running in + CI, which is what stops this recurring; it is turned on above. + +- **``Gamma`` no longer offers probability plotting as a fit method.** + ``Gamma.fit(x, how="MPP")`` now raises, joining ``Beta`` and + ``ExpoWeibull``, which already declined for the same reason. + + A probability plot works by rearranging the survival function so some + transform of the data falls on a straight line. For a Weibull, + ``log(-log S) = beta log x - beta log alpha`` — the axes do not depend + on the answer, so you can draw them before knowing anything. The + Gamma has no such rearrangement: its CDF is the regularised incomplete + gamma function and the shape sits *inside* that special function + rather than outside as an exponent. The only straight-line y-axis is + the inverse incomplete gamma, which needs the shape. To draw the axis + you need the answer; to get the answer you need the axis. + + The code broke the circle by guessing the shape from moments, drawing + the plot on that guess, and regressing. When the guess was off, the + axis was the wrong axis, the points were no longer straight on it, and + the regression fitted a line through a curve — returning a confident + wrong estimate rather than an error. An offset made it worse: the + shift distorts the low-``x`` end hardest, which is exactly where the + shape information lives. + + ``plot()`` is unaffected. It transforms with the *fitted* parameters, + so by the time the plot is drawn the axis is the right one — the + probability plot of an MLE-fitted Gamma remains a valid diagnostic. + Fitting is unchanged for MLE (the default), MSE and MOM. + + The 118-line ``Gamma.mpp`` override is deleted with it, which removes + the ``rr="x"`` mis-inversion and the censored-data ``LinAlgError`` from + #257 by making both paths unreachable. The MPP sweep in ``test_fit.py`` + now gates on each distribution's ``supports_mpp`` flag instead of a + hardcoded exclusion list, so it stays correct without editing. + +- **Nine wrong examples in the distribution docstrings.** Running + ``pytest --doctest-modules`` over ``distributions/`` gives 39 failures. + Thirty are the numpy-2 scalar repr (``np.float64(2.719...)`` against a + recorded ``2.719...``) and are cosmetic. Nine were not. + + Five documented outputs were simply wrong. ``Uniform.ff``'s example + called ``Uniform.sf``, and ``ExpoWeibull.cs``'s called + ``ExpoWeibull.sf`` -- in both the printed values were right for the + function being documented and wrong for the one being called, so the + example read as if the two were the same. ``LogLogistic.sf`` carried + values from some other parameterisation entirely (0.622 where the + answer is 0.988), ``LogLogistic.mean(3, 4)`` claimed ``3`` against + ``3.3322`` (the closed form is ``alpha (pi/beta) / sin(pi/beta)``), and + ``Exponential.qf`` had stale digits. + + The other four were the ``CustomDistribution`` example -- the Gompertz + walkthrough -- whose multi-line ``def`` used ``>>>`` where doctest + needs ``...``, so pasting it raised ``IndentationError``. + + In every case the code was right and the documentation was wrong, which + is the reassuring direction, but a reader checking their understanding + against these would have been misled. They accumulated precisely + because the doctests were not run, which is addressed above. + +- **Documented why a Turnbull fit does not equal a Kaplan-Meier fit.** + ``Turnbull.fit`` defaults to ``turnbull_estimator="Fleming-Harrington"`` + while ``KaplanMeier.fit`` is, unsurprisingly, KM. The Turnbull EM + recovers the same ``r`` and ``d`` either way; the three estimator + options then differ in how those become a survival curve. Comparing the + default against ``KaplanMeier`` and reading the gap as a defect is an + easy mistake — it is the one #260 was filed on, and the one made again + while checking whether #260 was still open. + + On ``x=[2,3,3,4,5,6], tl=[0,0,1,1,2,2]`` the survival at 2 is 0.750 + under KM, 0.765 under FH and 0.779 under NA. With the estimator matched, + Turnbull agrees with ``KaplanMeier`` to around 1e-9 on both ``sf`` and + ``cb``, across right-censored and left-truncated data. + + Only the KM option is the non-parametric MLE. Maximising the truncated + likelihood directly over the mass vector gives 0.750; FH's 0.765 scores + worse on that same likelihood, as an ``exp(-H)`` construction should. + FH is the default because it behaves better in the far tails and on + zero-inflated data (v0.8.0), not because it maximises anything. The + docstring now says all of this, and a test pins the three figures and + the NPMLE identity against a brute-force maximisation. + + No behaviour change. + v0.19.0 (4 August 2026) ----------------------- diff --git a/docs/requirements.txt b/docs/requirements.txt deleted file mode 100644 index 22d611ee..00000000 --- a/docs/requirements.txt +++ /dev/null @@ -1,17 +0,0 @@ -# Documentation build toolchain (used by ReadTheDocs and local builds). -# -# This is a coherent, fully pinned set validated by a complete docs -# build (every jupyter-execute cell runs) in a clean virtualenv on top -# of `pip install .`. Bump versions together and re-validate with: -# python -m sphinx -b html docs docs/_build/html -# -# ipykernel is capped below 7: jupyter-sphinx notebook execution hangs -# or crashes (zmq "Socket operation on non-socket") against the -# ipykernel 7 line. -matplotlib>=3.10 -sphinx==8.2.3 -sphinx_rtd_theme==3.1.0 -sphinx-copybutton==0.5.2 -sphinx-notfound-page==1.1.0 -jupyter-sphinx==0.5.3 -ipykernel==6.31.0 diff --git a/pyproject.toml b/pyproject.toml index 98f6d8ae..67545eea 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -39,6 +39,29 @@ tests = [ # (MongoDB compatibility); pymongo provides the bson package. "pymongo>=4", ] +docs = [ + # The documentation build toolchain, used by Read the Docs and by a + # local ``sphinx-build``. Fully pinned rather than bounded: the docs + # execute every ``.. jupyter-execute::`` cell during the build, so a + # toolchain change can break the build in ways a version range would + # let through silently. Bump these together and re-validate with a + # complete build: + # + # pip install -e ".[docs]" + # python -m sphinx -b html docs docs/_build/html + # + # matplotlib is not repeated here -- it is a runtime dependency of + # the package itself, which the docs install alongside these. + "sphinx==8.2.3", + "sphinx_rtd_theme==3.1.0", + "sphinx-copybutton==0.5.2", + "sphinx-notfound-page==1.1.0", + "jupyter-sphinx==0.5.3", + # Capped below 7: jupyter-sphinx notebook execution hangs or dies + # with a zmq "Socket operation on non-socket" against the ipykernel + # 7 line. + "ipykernel==6.31.0", +] [project.urls] Homepage = "https://github.com/derrynknife/SurPyval" @@ -49,6 +72,22 @@ where = ["."] [tool.setuptools.package-data] "*" = ["datasets/*.csv", "py.typed"] +# pytest +[tool.pytest.ini_options] +# Applied to the ``--doctest-modules`` run (see .github/workflows/actions.yml), +# which executes the ``>>>`` examples in the docstrings and compares their +# printed output. +# +# numpy chooses its own line breaks and column padding when it reprs an +# array, and both change with the width of the widest element. Without +# this an example is only correct for the exact terminal wrapping it was +# captured at, and a docstring wrapped to 79 columns by hand -- as most +# of these are -- can never match. +# +# The other half of the comparison is in conftest.py, which compares the +# numbers in the output as numbers rather than as text. +doctest_optionflags = "NORMALIZE_WHITESPACE" + # For Black pre-commit hook [tool.black] line-length = 79 diff --git a/scripts/check_all_pythons.py b/scripts/check_all_pythons.py new file mode 100644 index 00000000..5f052f14 --- /dev/null +++ b/scripts/check_all_pythons.py @@ -0,0 +1,195 @@ +#!/usr/bin/env python3 +"""Run the deployment checks locally, across every supported Python. + +Continuous integration runs the test suite only on the release pull +request into ``master`` (see docs/Contributing.rst). That keeps the +edit-review loop fast, at the cost of finding an interpreter-specific +failure at release time rather than when it lands. This script is the +other half of that trade: it runs what CI would have run, on all three +interpreters, before you push. + +It exists because the failure it guards against is not hypothetical. +The doctest numeric comparison added in 0.19.1 passed on 3.11 -- the +interpreter it was written on -- and failed on 3.12 and 3.13, because an +optimiser landed on a different last digit. Nothing short of actually +running the other interpreters would have found it. + + python scripts/check_all_pythons.py # 3.11, 3.12, 3.13 + python scripts/check_all_pythons.py 3.12 # just one + python scripts/check_all_pythons.py --skip-install # reuse as-is + +Environments live in ``.venvs/py3.X`` (git-ignored) and are reused +between runs, so only the first is slow. ``uv`` is used when it is +installed, because it makes the install step take seconds; otherwise +this falls back to ``venv`` and ``pip``, which works the same and takes +longer. + +Exits non-zero if any check on any interpreter fails. +""" + +from __future__ import annotations + +import argparse +import shutil +import subprocess +import sys +import time +from pathlib import Path + +ROOT = Path(__file__).resolve().parent.parent +VENVS = ROOT / ".venvs" +DEFAULT_VERSIONS = ["3.11", "3.12", "3.13"] + +# What CI runs, minus lint (which runs on every pull request anyway, so +# it is not what this script is for). Kept in the same order and with +# the same arguments as .github/workflows/actions.yml -- if they drift, +# this stops being a preview of CI and becomes its own thing. +CHECKS: list[tuple[str, list[str]]] = [ + ( + "test suite", + [ + "-m", + "pytest", + "-n", + "auto", + "-q", + "--ignore=surpyval/tests/alpha", + "--run-ml", + ], + ), + ( + "doctests", + [ + "-m", + "pytest", + "--doctest-modules", + "surpyval", + "-q", + "--ignore=surpyval/tests", + "--ignore=surpyval/alpha", + ], + ), + ( + "doctests (numeric forced)", + [ + "-m", + "pytest", + "--doctest-modules", + "surpyval", + "-q", + "--ignore=surpyval/tests", + "--ignore=surpyval/alpha", + "--doctest-force-numeric", + ], + ), +] + + +def run(cmd: list[str], **kwargs) -> int: + """Run a command, streaming its output, and return its exit status.""" + return subprocess.call(cmd, cwd=ROOT, **kwargs) + + +def interpreter_for(venv: Path) -> Path: + bin_dir = "Scripts" if sys.platform == "win32" else "bin" + exe = "python.exe" if sys.platform == "win32" else "python" + return venv / bin_dir / exe + + +def ensure_env(version: str, skip_install: bool) -> Path | None: + """Create (or reuse) the environment for ``version``. + + Returns the interpreter path, or None if the interpreter is not + available on this machine -- which is a thing to report, not a + thing to crash on. + """ + venv = VENVS / f"py{version}" + python = interpreter_for(venv) + uv = shutil.which("uv") + + if not python.exists(): + print(f"\n[{version}] creating {venv.relative_to(ROOT)}") + if uv: + created = run([uv, "venv", "--python", version, str(venv)]) + else: + base = shutil.which(f"python{version}") + if base is None: + print(f"[{version}] no python{version} on PATH -- skipping") + return None + created = run([base, "-m", "venv", str(venv)]) + if created != 0 or not python.exists(): + print(f"[{version}] could not create an environment -- skipping") + return None + + if not skip_install: + print(f"[{version}] installing") + # pytest-xdist is not in the tests extra; it is a development + # tool, and the suite's ``-n auto`` needs it. + target = ["-e", ".[tests]", "pytest-xdist"] + if uv: + installed = run( + [uv, "pip", "install", "-q", "--python", str(python), *target] + ) + else: + installed = run( + [str(python), "-m", "pip", "install", "-q", *target] + ) + if installed != 0: + print(f"[{version}] install failed -- skipping") + return None + + return python + + +def main() -> int: + parser = argparse.ArgumentParser( + description="Run the CI checks locally on every supported Python." + ) + parser.add_argument( + "versions", + nargs="*", + default=DEFAULT_VERSIONS, + help="Python versions to check (default: %s)" + % " ".join(DEFAULT_VERSIONS), + ) + parser.add_argument( + "--skip-install", + action="store_true", + help="reuse the environments as they are, without reinstalling", + ) + args = parser.parse_args() + + VENVS.mkdir(exist_ok=True) + results: list[tuple[str, str, str, float]] = [] + + for version in args.versions: + python = ensure_env(version, args.skip_install) + if python is None: + results.append((version, "environment", "UNAVAILABLE", 0.0)) + continue + for name, cmd in CHECKS: + print(f"\n[{version}] {name}") + started = time.monotonic() + status = run([str(python), *cmd]) + elapsed = time.monotonic() - started + results.append( + (version, name, "ok" if status == 0 else "FAILED", elapsed) + ) + + print("\n" + "=" * 60) + failed = False + for version, name, status, elapsed in results: + if status != "ok": + failed = True + print(f"{version:6s} {name:28s} {status:12s} {elapsed:6.1f}s") + print("=" * 60) + + if failed: + print("\nSomething failed. Do not push.") + return 1 + print("\nAll checks passed on every interpreter.") + return 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/surpyval/degradation/degradation_analysis.py b/surpyval/degradation/degradation_analysis.py index 2652d4bf..0197fd5c 100644 --- a/surpyval/degradation/degradation_analysis.py +++ b/surpyval/degradation/degradation_analysis.py @@ -1193,8 +1193,8 @@ class DegradationAnalysis_: Censored Units : 0 Life Distribution : Weibull Parameters : - alpha: 441.4780882117898 - beta: 6.987078993008337 + alpha: 441.47809611105606 + beta: 6.987078889297555 >>> model.pseudo_failure_times array([451.61290323, 500. , 318.18181818, 378.37837838]) """ diff --git a/surpyval/recurrent/parametric/cox_lewis.py b/surpyval/recurrent/parametric/cox_lewis.py index f1b3270a..e66f9cdf 100644 --- a/surpyval/recurrent/parametric/cox_lewis.py +++ b/surpyval/recurrent/parametric/cox_lewis.py @@ -27,19 +27,19 @@ class CoxLewis(NHPPFitter): Process : Cox-Lewis Fitted by : MLE Parameters : - alpha: 0.3848528712360503 - beta: 0.1939447728437042 + alpha: 0.384812737762836 + beta: 0.19396672109211047 >>> model.cif([1, 2, 3, 4, 5, 6]) - array([ 1.6215655 , 3.59019342, 5.98016527, 8.88166096, 12.40416155, - 16.68058024]) + array([ 1.62151879, 3.59013322, 5.98014113, 8.88174429, 12.40445268, + 16.6812175 ]) >>> >>> model.iif([1, 2, 3, 4, 5, 6]) - array([1.78389227, 2.16569736, 2.62921991, 3.19194983, 3.87512041, - 4.70450946]) + array([1.78385983, 2.16570551, 2.62928751, 3.19210196, 3.87539016, + 4.70494021]) >>> >>> model.inv_cif([1, 2, 3, 4, 5, 6]) - array([0.63923607, 1.20789182, 1.72001505, 2.18583659, 2.61303902, - 3.00753845]) + array([0.63925589, 1.20792021, 1.72004461, 2.18586234, 2.61305756, + 3.00754742]) """ def __init__(self): diff --git a/surpyval/recurrent/parametric/crow_amsaa.py b/surpyval/recurrent/parametric/crow_amsaa.py index ef66cf05..2bd8eac7 100644 --- a/surpyval/recurrent/parametric/crow_amsaa.py +++ b/surpyval/recurrent/parametric/crow_amsaa.py @@ -27,19 +27,19 @@ class CrowAMSAA(NHPPFitter): Process : Crow-AMSAA Fitted by : MLE Parameters : - alpha: 3129.2801848331596 - beta: 1.239258986741094 + alpha: 913.8466210685444 + beta: 1.4781707110680866 >>> model.cif([1, 2, 3, 4, 5, 6]) - array([4.65842366e-05, 1.09974784e-04, 1.81767316e-04, 2.59625443e-04, - 3.42329130e-04, 4.29111278e-04]) + array([4.20072057e-05, 1.17030084e-04, 2.13103439e-04, 3.26040266e-04, + 4.53440995e-04, 5.93696079e-04]) >>> >>> model.iif([1, 2, 3, 4, 5, 6]) - array([5.77299348e-05, 6.81436206e-05, 7.50855947e-05, 8.04357921e-05, - 8.48468915e-05, 8.86300027e-05]) + array([6.20938211e-05, 8.64952211e-05, 1.05001087e-04, 1.20485793e-04, + 1.34052640e-04, 1.46264026e-04]) >>> >>> model.inv_cif([1, 2, 3, 4, 5, 6]) - array([ 3129.2801836 , 5474.64210552, 7593.6351178 , 9577.82762328, - 11467.45338053, 13284.98316248]) + array([ 913.84662107, 1460.57434899, 1921.54912724, 2334.39329941, + 2714.78099355, 3071.15581638]) """ def __init__(self): diff --git a/surpyval/recurrent/parametric/duane.py b/surpyval/recurrent/parametric/duane.py index 71f7963f..0974d1ca 100644 --- a/surpyval/recurrent/parametric/duane.py +++ b/surpyval/recurrent/parametric/duane.py @@ -27,19 +27,19 @@ class Duane(NHPPFitter): Process : Duane Fitted by : MLE Parameters : - alpha: 1.2392945732132952 - b: 4.6568641229556424e-05 + alpha: 1.478202089169939 + b: 4.199455086392048e-05 >>> model.cif([1, 2, 3, 4, 5, 6]) - array([4.65686412e-05, 1.09940677e-04, 1.81713565e-04, 2.59551323e-04, - 3.42234115e-04, 4.28994958e-04]) + array([4.19945509e-05, 1.16997373e-04, 2.13046585e-04, 3.25956224e-04, + 4.53327287e-04, 5.93550595e-04]) >>> >>> model.iif([1, 2, 3, 4, 5, 6]) - array([5.77122644e-05, 6.81244421e-05, 7.50655449e-05, 8.04151364e-05, - 8.48257764e-05, 8.86085206e-05]) + array([6.20764328e-05, 8.64728804e-05, 1.04975302e-04, 1.20457293e-04, + 1.34021869e-04, 1.46231288e-04]) >>> >>> model.inv_cif([1, 2, 3, 4, 5, 6]) - array([ 3129.4028404 , 5474.76881015, 7593.73955973, 9577.89554535, - 11467.47544342, 13284.95262974]) + array([ 913.90063856, 1460.64614435, 1921.63239303, 2334.48481047, + 2714.87871658, 3071.25832646]) """ def __init__(self): diff --git a/surpyval/recurrent/parametric/hpp.py b/surpyval/recurrent/parametric/hpp.py index 4bbc1ea9..43cbdf62 100644 --- a/surpyval/recurrent/parametric/hpp.py +++ b/surpyval/recurrent/parametric/hpp.py @@ -34,18 +34,17 @@ class HPP(CountingProcess): Process : Homogeneous Poisson Process Fitted by : MLE Parameters : - lambda: 0.000498450145693719 + lambda: 0.0023047023327236213 >>> model.cif([1, 2, 3, 4, 5, 6]) - array([0.00049845, 0.0009969 , 0.00149535, 0.0019938 , 0.00249225, - 0.0029907 ]) + array([0.0023047 , 0.0046094 , 0.00691411, 0.00921881, 0.01152351, + 0.01382821]) >>> >>> model.iif([1, 2, 3, 4, 5, 6]) - array([0.00049845, 0.00049845, 0.00049845, 0.00049845, 0.00049845, - 0.00049845]) + array([0.0023047, 0.0023047, 0.0023047, 0.0023047, 0.0023047, 0.0023047]) >>> >>> model.inv_cif([1, 2, 3, 4, 5, 6]) - array([ 2006.21869336, 4012.43738672, 6018.65608009, 8024.87477345, - 10031.09346681, 12037.31216017]) + array([ 433.89551258, 867.79102516, 1301.68653774, 1735.58205032, + 2169.4775629 , 2603.37307548]) """ def __init__(self): diff --git a/surpyval/recurrent/regression/proportional_intensity.py b/surpyval/recurrent/regression/proportional_intensity.py index 3844f2fa..4a4d4035 100644 --- a/surpyval/recurrent/regression/proportional_intensity.py +++ b/surpyval/recurrent/regression/proportional_intensity.py @@ -37,10 +37,10 @@ class ProportionalIntensityModel( >>> c = data['arrest'].values >>> Z = data[["fin", "age", "race", "wexp", "mar", "paro", "prio"]].values >>> model = ProportionalIntensityNHPP.fit(x, Z, c, dist=CrowAMSAA) - >>> type(model) - surpyval.recurrent.regression.proportional_intensity.ProportionalIntensityModel + >>> type(model).__name__ + 'ProportionalIntensityModel' >>> model.cif([1, 2, 3], Z.mean(axis=0)) - array([0.00625402, 0.04304137, 0.13302238]) + array([8.84210972e-07, 2.79074784e-05, 2.10220821e-04]) """ def __repr__(self): diff --git a/surpyval/recurrent/renewal/generalized_one_renewal.py b/surpyval/recurrent/renewal/generalized_one_renewal.py index b54e148d..a9d1937c 100644 --- a/surpyval/recurrent/renewal/generalized_one_renewal.py +++ b/surpyval/recurrent/renewal/generalized_one_renewal.py @@ -60,16 +60,16 @@ class GeneralizedOneRenewal(RenewalFitMixin): ========================= Distribution : Weibull Fitted by : MLE - Restoration Factor : -0.1730179893443181 + Restoration Factor : -0.1730184624683848 Parameters : - alpha: 1.3919016662855024 - beta: 5.008872636271443 + alpha: 1.3919045968817332 + beta: 5.008861189641614 >>> >>> np.random.seed(0) >>> np_model = model.count_terminated_simulation(len(x), 5000) >>> np_model.mcf(np.array([1, 2, 3, 4, 5, 6])) - array([0.1696 , 1.181 , 2.287 , 3.6696 , 5.58237921, - 8.54474127]) + array([0.1696 , 1.181 , 2.287 , 3.6694 , 5.58237925, + 8.54474531]) """ @staticmethod @@ -208,10 +208,10 @@ def fit_from_recurrent_data(self, data, dist=Weibull, init=None): ========================= Distribution : Weibull Fitted by : MLE - Restoration Factor : 0.4270960618530103 + Restoration Factor : 0.3402789091696592 Parameters : - alpha: 1.3494830373118245 - beta: 2.7838386997223212 + alpha: 1.4115217370254167 + beta: 3.5499343659245564 """ self._check_dist_eligible(dist) validate_renewal_censoring(data.c, type(self).__name__) @@ -290,10 +290,10 @@ def fit(self, x, i=None, c=None, n=None, dist=Weibull, init=None): ========================= Distribution : Weibull Fitted by : MLE - Restoration Factor : 0.4270960618530103 + Restoration Factor : 0.3402789091696592 Parameters : - alpha: 1.3494830373118245 - beta: 2.7838386997223212 + alpha: 1.4115217370254167 + beta: 3.5499343659245564 """ data = handle_xicn(x, i, c, n) return self.fit_from_recurrent_data(data, dist=dist, init=init) @@ -325,10 +325,10 @@ def fit_from_parameters(self, params, q, dist=Weibull): >>> from surpyval.recurrent import GeneralizedOneRenewal >>> >>> model = GeneralizedOneRenewal.fit_from_parameters( - [10, 2], - 0.2, - dist=Weibull - ) + ... [10, 2], + ... 0.2, + ... dist=Weibull + ... ) >>> model G1 Renewal SurPyval Model ========================= diff --git a/surpyval/recurrent/renewal/generalized_renewal.py b/surpyval/recurrent/renewal/generalized_renewal.py index 71210fa7..eaa40bcf 100644 --- a/surpyval/recurrent/renewal/generalized_renewal.py +++ b/surpyval/recurrent/renewal/generalized_renewal.py @@ -60,10 +60,10 @@ class GeneralizedRenewal(RenewalFitMixin): Distribution : Weibull Fitted by : MLE Kijima Type : i - Restoration Factor : 0.1573211400037486 + Restoration Factor : 0.15732122999163628 Parameters : - alpha: 1.261338468404201 - beta: 8.93900788677076 + alpha: 1.261337933121844 + beta: 8.93902321971521 >>> >>> np.random.seed(0) >>> np_model = model.count_terminated_simulation(len(x), 5000) @@ -259,10 +259,10 @@ def fit_from_recurrent_data( Distribution : Weibull Fitted by : MLE Kijima Type : i - Restoration Factor : 1.594694243423234e-11 + Restoration Factor : 1.3316262291443964e-16 Parameters : - alpha: 2.399029078569064 - beta: 2.753920439616154 + alpha: 2.399029668688425 + beta: 2.753920042066547 """ validate_renewal_censoring(data.c, type(self).__name__) reject_left_truncation(data, type(self).__name__) @@ -344,10 +344,10 @@ def fit( Distribution : Weibull Fitted by : MLE Kijima Type : i - Restoration Factor : 1.594694243423234e-11 + Restoration Factor : 1.3316262291443964e-16 Parameters : - alpha: 2.399029078569064 - beta: 2.753920439616154 + alpha: 2.399029668688425 + beta: 2.753920042066547 """ data = handle_xicn(x, i, c, n) return self.fit_from_recurrent_data(data, dist, kijima, init=init) @@ -381,10 +381,10 @@ def fit_from_parameters(self, params, q, kijima="i", dist=Weibull): >>> from surpyval.recurrent import GeneralizedRenewal >>> >>> model = GeneralizedRenewal.fit_from_parameters( - [10, 2], - 0.2, - dist=Normal - ) + ... [10, 2], + ... 0.2, + ... dist=Normal + ... ) >>> model Generalized Renewal SurPyval Model ================================== diff --git a/surpyval/tests/test_doctest_checker.py b/surpyval/tests/test_doctest_checker.py new file mode 100644 index 00000000..9193ff3f --- /dev/null +++ b/surpyval/tests/test_doctest_checker.py @@ -0,0 +1,165 @@ +"""The numeric fallback used by the ``--doctest-modules`` run. + +``conftest._numerically_equal`` decides whether two blocks of doctest +output that differ as text are the same as numbers. It is the mechanism +that lets the docstring examples record real, untrimmed output while +still passing on every supported Python, so it needs its own tests: too +strict and the doctest step fails on a toolchain difference, too loose +and it stops catching the documentation drift it was added to catch. + +The "tolerated" cases below are not invented. Each is a real pair of +outputs seen for the same example on two different Pythons in CI. +""" + +import doctest + +import pytest + +from conftest import _numerically_equal + +# The ``ProportionalIntensityHPP`` example, verbatim: the documented +# block on the left, the block CI produced on the right. Kept whole +# rather than reduced to one coefficient because the thing that broke +# here was not a number -- it was the ```` markers, which the +# text comparison strips from the expected output before matching and +# which the numeric fallback has to strip too. A test on an extracted +# line would have gone on passing. +HPP_DOCUMENTED = """\ +Base Rate Parameters: + lambda : 0.012395105741757225 + +Covariate Coefficients: + beta_0 : 0.06397367067847898 + beta_1 : 0.011491178797116433 + beta_2 : -0.02147901865302258 + +""" + +HPP_OBSERVED = """\ +Base Rate Parameters: + lambda : 0.012395109943236718 + +Covariate Coefficients: + beta_0 : 0.06397361219411789 + beta_1 : 0.011491197469342556 + beta_2 : -0.021479544898597686 + +""" + +# (documented output, output seen on another Python) +TOLERATED = [ + pytest.param( + " alpha: 913.84662107\n beta: 1.4781707110680866\n", + " alpha: 913.8468395959314\n beta: 1.4781709287931744\n", + id="crow-amsaa-fit", + ), + pytest.param( + # The loosest genuine disagreement found: 9e-4 relative on ``b``. + " alpha: 1.478202089169939\n b: 4.199455086392048e-05\n", + " alpha: 1.4781024854527343\n b: 4.203204430005199e-05\n", + id="duane-fit", + ), + pytest.param( + " beta_2 : -0.02147901865302258\n", + " beta_2 : -0.021479544898597686\n", + id="hpp-proportional-intensity-coefficient", + ), + pytest.param( + # A restoration factor whose true value is zero; the optimiser + # stops on whatever mantissa it stops on. Only ``abs_tol`` can + # rescue this one. + "Restoration Factor : 1.3316262291443964e-16\n", + "Restoration Factor : 1.1551809284521243e-16\n", + id="numerical-zero", + ), + pytest.param( + "1.8227536487527594\n", + "1.822753648752769\n", + id="scipy-special-last-digits", + ), + pytest.param( + # numpy re-pads the columns when an element gets shorter, so the + # closing bracket moves. + " alpha: [ 6.32508961 17.37701969]\n", + " alpha: [ 6.3250866 17.377018 ]\n", + id="array-repadded", + ), + pytest.param(HPP_DOCUMENTED, HPP_OBSERVED, id="repr-with-blank-lines"), +] + +# Every one of these is a real defect this sweep found in the docstrings, +# or the class of defect it found. None may be forgiven. +REJECTED = [ + pytest.param( + "3\n", + "3.332162203618775\n", + id="value-from-another-parameterisation", + ), + pytest.param( + "11.229\n", "10.533288486847923\n", id="documented-variance-was-wrong" + ), + pytest.param( + "array([1, 2, 3, 4, 5])\n", + "array([1., 2., 3., 4., 5.])\n", + id="dtype-changed", + ), + pytest.param( + "array([0.83333333, 0.66666667])\n", + "array([0.16666667, 0.33333333])\n", + id="example-called-the-wrong-function", + ), + pytest.param( + "array([1.0, 2.0])\n", "array([1.0, 2.0, 3.0])\n", id="shape-changed" + ), + pytest.param("alpha: 1.0\n", "beta: 1.0\n", id="different-label"), + pytest.param("array([1.0])\n", "3.0\n", id="scalar-against-array"), + pytest.param("", "1.0\n", id="expected-nothing"), + pytest.param("1.0\n", "", id="produced-nothing"), + pytest.param( + # 2e-3 relative, twice the tolerance. + "8.929795115692489\n", + "8.947654505923874\n", + id="drifted-past-the-tolerance", + ), +] + + +@pytest.mark.parametrize("want, got", TOLERATED) +def test_the_same_number_under_a_different_toolchain_is_accepted(want, got): + assert _numerically_equal(want, got) + + +@pytest.mark.parametrize("want, got", REJECTED) +def test_a_real_difference_is_still_a_failure(want, got): + assert not _numerically_equal(want, got) + + +def test_output_with_no_numbers_is_not_silently_accepted(): + # Nothing to compare numerically, so the fallback must decline and + # leave the verdict to the ordinary text comparison. + assert not _numerically_equal("Fitted by : MLE\n", "Fitted by : MPS\n") + assert not _numerically_equal("Fitted by : MLE\n", "Fitted by : MLE\n") + + +def test_the_checker_is_installed_on_the_doctest_base_class(): + # The examples' precision depends on this patch being live for the + # whole doctest run, including under pytest's own subclass. + checker = doctest.OutputChecker() + assert checker.check_output("1.0000000\n", "1.0000001\n", 0) + assert not checker.check_output("1.0\n", "2.0\n", 0) + + +def test_the_whole_comparison_path_handles_blank_lines(): + # Not the fallback in isolation but the method doctest actually + # calls, with the flags the doctest step actually runs under. The + # ```` handling lives in the text comparison, so only + # this route proves the two halves agree about it. + checker = doctest.OutputChecker() + assert checker.check_output( + HPP_DOCUMENTED, HPP_OBSERVED, doctest.NORMALIZE_WHITESPACE + ) + assert not checker.check_output( + HPP_DOCUMENTED, + HPP_OBSERVED.replace("0.0123951", "0.0198765"), + doctest.NORMALIZE_WHITESPACE, + ) diff --git a/surpyval/tests/univariate/nonparametric/test_turnbull.py b/surpyval/tests/univariate/nonparametric/test_turnbull.py index e72cdedf..dfe487ef 100644 --- a/surpyval/tests/univariate/nonparametric/test_turnbull.py +++ b/surpyval/tests/univariate/nonparametric/test_turnbull.py @@ -367,3 +367,67 @@ def test_max_iter_zero_raises(): surpyval.Turnbull.fit( np.array([1.0, 2.0, 3.0]), c=np.array([0, 1, 0]), max_iter=0 ) + + +def test_only_the_km_option_is_the_npmle(): + # The three turnbull_estimator options are not three ways of computing + # one number: the EM recovers the same r and d, and they then differ in + # how those become a survival curve. Only Kaplan-Meier is the + # non-parametric MLE; Nelson-Aalen and Fleming-Harrington are exp(-H) + # constructions that are not maximising anything. + # + # This pins the figures quoted in the ``turnbull_estimator`` docstring, + # which exist because comparing a default (FH) Turnbull fit against + # KaplanMeier and reading the gap as a defect is an easy mistake -- it + # is the mistake #260 was originally filed on. + from scipy.optimize import minimize + + x = np.array([2.0, 3.0, 3.0, 4.0, 5.0, 6.0]) + tl = np.array([0.0, 0.0, 1.0, 1.0, 2.0, 2.0]) + support = np.array([2.0, 3.0, 4.0, 5.0, 6.0]) + + def neg_ll(u): + # masses on the support, softmax-parameterised so they stay simplex + p = np.exp(u - u.max()) + p = p / p.sum() + total = 0.0 + for xi, ti in zip(x, tl): + mass = p[support == xi].sum() + at_risk = p[support > ti].sum() # (entry, exit] + if mass <= 0 or at_risk <= 0: + return 1e6 + total += np.log(mass) - np.log(at_risk) + return -total + + best = None + for seed in range(6): + res = minimize( + neg_ll, + np.random.default_rng(seed).normal(size=support.size), + method="Nelder-Mead", + options={"maxiter": 40000, "fatol": 1e-14, "xatol": 1e-12}, + ) + if best is None or res.fun < best.fun: + best = res + + p = np.exp(best.x - best.x.max()) + p = p / p.sum() + npmle_sf2 = p[support > 2.0].sum() + assert npmle_sf2 == pytest.approx(0.75, abs=1e-5) + + got = { + est: float( + np.ravel( + surpyval.Turnbull.fit(x, tl=tl, turnbull_estimator=est).sf(2.0) + )[0] + ) + for est in ("Kaplan-Meier", "Fleming-Harrington", "Nelson-Aalen") + } + + # KM reaches the NPMLE; the other two are elsewhere, and are ordered. + assert got["Kaplan-Meier"] == pytest.approx(npmle_sf2, abs=1e-5) + assert got["Fleming-Harrington"] == pytest.approx(0.7652, abs=1e-3) + assert got["Nelson-Aalen"] == pytest.approx(0.7788, abs=1e-3) + assert ( + got["Kaplan-Meier"] < got["Fleming-Harrington"] < got["Nelson-Aalen"] + ) diff --git a/surpyval/tests/univariate/parametric/test_distribution_fixes.py b/surpyval/tests/univariate/parametric/test_distribution_fixes.py index 54580272..b8ffe64c 100644 --- a/surpyval/tests/univariate/parametric/test_distribution_fixes.py +++ b/surpyval/tests/univariate/parametric/test_distribution_fixes.py @@ -44,18 +44,31 @@ def test_exponential_offset_mpp_rr_x_inversion(): def test_gamma_censored_mpp_rr_x_does_not_crash(): + # Was a LinAlgError: rr="x" regressed the filtered y against the raw + # x whenever censoring filtered any point (#257). Gamma declines MPP + # entirely now (#158), so the crashing path is unreachable -- the + # refusal arrives before any regression is attempted. np.random.seed(1) x = Gamma.random(300, 3, 2) c = (x > 2.5).astype(int) - m = Gamma.fit(np.minimum(x, 2.5), c=c, how="MPP", rr="x") - assert np.all(np.isfinite(m.params)) + with pytest.raises(ValueError, match="does not work with probability"): + Gamma.fit(np.minimum(x, 2.5), c=c, how="MPP", rr="x") @pytest.mark.parametrize("rr", ["x", "y"]) def test_gamma_offset_mpp_recovers_offset(rr): + # Gamma no longer offers MPP at all (#158): the probability plot's + # own y-axis is the inverse incomplete gamma, which needs the shape + # being estimated, so the fit regressed against an axis built from a + # guess and returned a confident wrong answer. It now refuses. np.random.seed(2) x = Gamma.random(300, 3, 2) + 10 - m = Gamma.fit(x, how="MPP", rr=rr, offset=True) + with pytest.raises(ValueError, match="does not work with probability"): + Gamma.fit(x, how="MPP", rr=rr, offset=True) + + # The offset recovery this test existed to protect still holds under + # the estimators Gamma does support. + m = Gamma.fit(x, offset=True) assert m.gamma == pytest.approx(10.0, abs=1.5) assert m.params[0] == pytest.approx(3.0, rel=0.5) diff --git a/surpyval/tests/univariate/parametric/test_fit.py b/surpyval/tests/univariate/parametric/test_fit.py index cc7e5706..b2d63f20 100644 --- a/surpyval/tests/univariate/parametric/test_fit.py +++ b/surpyval/tests/univariate/parametric/test_fit.py @@ -246,7 +246,10 @@ def test_mle_convergence_small(dist, random_parameters, kind): "dist,random_parameters,rr", generate_mpp_test_cases(), ids=idfunc ) def test_mpp(dist, random_parameters, rr): - if dist not in [Beta, ExpoWeibull]: + # Gate on the distribution's own flag rather than a hardcoded + # list: Beta, ExpoWeibull and Gamma all lack a linearising + # probability plot, and a new one should not need this edited. + if dist.supports_mpp: for n in FIT_SIZES: test_params = [] tol = 0.025 diff --git a/surpyval/tests/univariate/parametric/test_offset_divergence.py b/surpyval/tests/univariate/parametric/test_offset_divergence.py index fe14c938..18284712 100644 --- a/surpyval/tests/univariate/parametric/test_offset_divergence.py +++ b/surpyval/tests/univariate/parametric/test_offset_divergence.py @@ -13,6 +13,7 @@ """ import numpy as np +import pytest from scipy.stats import wasserstein_distance from surpyval import Gamma, Rayleigh @@ -64,19 +65,29 @@ def _summary(dist, true_params, how, seed=0): } -def test_mpp_offset_gamma_parameters_recovered(): +def test_mpp_offset_gamma_is_refused(): """MPP offset on Gamma used to land on an absurd parameter tuple (huge shape, gamma far off) that merely mimicked the true distribution. This was the divergence this module existed to pin - down; #257 fixed the initialiser (multi-started shape search) and - the rr="x" inversion, so the *parameters* are now recovered too.""" - s = _summary(Gamma, (3.0, 2.0), how="MPP") + down. #257 fixed the initialiser and the rr="x" inversion, which + recovered the parameters, but the underlying circularity remained: + the plot's y-axis is the inverse incomplete gamma, so it has to be + drawn from a guess at the very shape being estimated. Gamma now + declines MPP outright (#158) rather than returning a confident + answer off a wrong axis.""" + x = Gamma.random(500, 3.0, 2.0) + TRUE_GAMMA + with pytest.raises(ValueError, match="does not work with probability"): + Gamma.fit(x, how="MPP", offset=True) + + +def test_mle_offset_gamma_parameters_recovered(): + """The recovery the MPP test above was really protecting, under an + estimator Gamma supports.""" + s = _summary(Gamma, (3.0, 2.0), how="MLE") - # The parameters are now close to the truth. assert abs(s["fit"].gamma - TRUE_GAMMA) < 1.0, s["fit"].gamma assert s["max_param_rel_err"] < 0.20 # within 20% - # ...and the distribution remains essentially exact. assert s["mean_rel_err"] < 0.01 # mean within 1% assert s["median_rel_err"] < 0.03 # median within 3% assert s["std_rel_err"] < 0.10 # spread within 10% diff --git a/surpyval/univariate/information_criteria.py b/surpyval/univariate/information_criteria.py index 7db588f8..d45960c3 100644 --- a/surpyval/univariate/information_criteria.py +++ b/surpyval/univariate/information_criteria.py @@ -59,7 +59,7 @@ def neg_ll(self) -> float: >>> x = Weibull.random(100, 10, 3) >>> model = Weibull.fit(x) >>> model.neg_ll() - 262.52685642385734 + 262.52685642390634 """ if getattr(self, "data", None) is None: raise ValueError("Must have been fit with data") @@ -88,7 +88,7 @@ def bic(self) -> float: >>> x = Weibull.random(100, 10, 3) >>> model = Weibull.fit(x) >>> model.bic() - 534.2640532196908 + np.float64(534.2640532197888) References ---------- @@ -124,7 +124,7 @@ def aic(self) -> float: >>> x = Weibull.random(100, 10, 3) >>> model = Weibull.fit(x) >>> model.aic() - 529.0537128477147 + 529.0537128478127 """ if hasattr(self, "_aic"): return self._aic @@ -152,7 +152,7 @@ def aic_c(self) -> float: >>> x = Weibull.random(100, 10, 3) >>> model = Weibull.fit(x) >>> model.aic_c() - 529.1774241879209 + np.float64(529.1774241880189) """ if hasattr(self, "_aic_c"): return self._aic_c diff --git a/surpyval/univariate/nonparametric/nonparametric.py b/surpyval/univariate/nonparametric/nonparametric.py index b9f41ad1..fadbeed5 100755 --- a/surpyval/univariate/nonparametric/nonparametric.py +++ b/surpyval/univariate/nonparametric/nonparametric.py @@ -686,7 +686,7 @@ def mean(self, tau: float | None = None) -> float: >>> x = np.array([1, 2, 3, 4, 5]) >>> model = KaplanMeier.fit(x) >>> model.mean() - 3.0 + 3.0000000000000004 """ if np.min(self.x) < 0: raise ValueError( diff --git a/surpyval/univariate/nonparametric/nonparametric_fitter.py b/surpyval/univariate/nonparametric/nonparametric_fitter.py index ecec6d77..e9011f8f 100755 --- a/surpyval/univariate/nonparametric/nonparametric_fitter.py +++ b/surpyval/univariate/nonparametric/nonparametric_fitter.py @@ -117,6 +117,26 @@ def fit( KM, NA, or FH estimator with the Turnbull estimates of r, and d. Defaults to FH. + **This default is why a Turnbull fit does not equal a + KaplanMeier fit on data both can handle.** The Turnbull EM + recovers the same ``r`` and ``d``; the three options then differ + in how they turn those into a survival curve, so the difference + is the estimator, not the data or the EM. On + ``x=[2,3,3,4,5,6], tl=[0,0,1,1,2,2]`` the survival at 2 is + 0.750 under KM, 0.765 under FH and 0.779 under NA. Pass + ``turnbull_estimator='Kaplan-Meier'`` to compare like with like + -- it then agrees with :code:`KaplanMeier` to around 1e-9 on + both ``sf`` and ``cb``, on right-censored and left-truncated + data alike. + + Only the KM option is the non-parametric MLE. NA and FH are + ``exp(-H)`` constructions and are not trying to maximise the + likelihood: on the data above the direct NPMLE of the truncated + likelihood is 0.750, and FH's 0.765 scores worse on it by + design. FH is the default because it behaves better than KM in + the far tails and on zero-inflated data (see the v0.8.0 notes), + not because it is the maximum likelihood answer. + tol : float, optional Turnbull only. The EM stops once the largest change in any interval's probability mass falls below this. Defaults to 1e-10. diff --git a/surpyval/univariate/nonparametric/plotting_positions.py b/surpyval/univariate/nonparametric/plotting_positions.py index 254ac29f..976ea04e 100755 --- a/surpyval/univariate/nonparametric/plotting_positions.py +++ b/surpyval/univariate/nonparametric/plotting_positions.py @@ -97,7 +97,7 @@ def plotting_positions( Examples -------- - >>> from surpyval.nonparametric import plotting_positions + >>> from surpyval.univariate.nonparametric import plotting_positions >>> import numpy as np >>> x = np.array([1, 2, 3, 4, 5, 6, 7, 8]) >>> x, r, d, F = plotting_positions(x, heuristic="Filliben") diff --git a/surpyval/univariate/nonparametric/success_run.py b/surpyval/univariate/nonparametric/success_run.py index a115b70c..5d27b6f3 100755 --- a/surpyval/univariate/nonparametric/success_run.py +++ b/surpyval/univariate/nonparametric/success_run.py @@ -32,7 +32,7 @@ def success_run(n, confidence=None, alpha=None): >>> from surpyval import success_run >>> success_run(10) - 0.7411344491069477 + np.float64(0.7411344491069477) """ if confidence and alpha: raise ValueError("Only one of confidence or alpha can be specified") diff --git a/surpyval/univariate/parametric/distributions/beta.py b/surpyval/univariate/parametric/distributions/beta.py index 670cf66e..b0ccca59 100755 --- a/surpyval/univariate/parametric/distributions/beta.py +++ b/surpyval/univariate/parametric/distributions/beta.py @@ -333,7 +333,7 @@ def moment(self, n, alpha, beta): -------- >>> from surpyval import Beta >>> Beta.moment(2, 3, 4) - 0.21428571428571427 + np.float64(0.2142857142857143) """ return np.exp(abetaln(n + alpha, beta) - abetaln(alpha, beta)) @@ -369,7 +369,7 @@ def entropy(self, alpha, beta): -------- >>> from surpyval import Beta >>> Beta.entropy(3, 4) - -0.3443445622221013 + np.float64(-0.3443445622221013) """ return ( abetaln(alpha, beta) diff --git a/surpyval/univariate/parametric/distributions/beta4.py b/surpyval/univariate/parametric/distributions/beta4.py index 815d5a66..3b2a6f89 100644 --- a/surpyval/univariate/parametric/distributions/beta4.py +++ b/surpyval/univariate/parametric/distributions/beta4.py @@ -411,7 +411,7 @@ def moment(self, m, alpha, beta, a, b): -------- >>> from surpyval import Beta4 >>> Beta4.moment(1, 3, 4, 2, 3) - 2.4285714285714284 + np.float64(2.428571428571429) """ scale = b - a total = 0.0 diff --git a/surpyval/univariate/parametric/distributions/binomial.py b/surpyval/univariate/parametric/distributions/binomial.py index 5038e049..ee6b8532 100644 --- a/surpyval/univariate/parametric/distributions/binomial.py +++ b/surpyval/univariate/parametric/distributions/binomial.py @@ -64,7 +64,7 @@ def df(self, x, n, p): -------- >>> from surpyval import Binomial >>> Binomial.df(2, 5, 0.3) - 0.3086999999999998 + np.float64(0.3086999999999998) """ return binom.pmf(x, n, p) @@ -97,7 +97,7 @@ def ff(self, x, n, p): -------- >>> from surpyval import Binomial >>> Binomial.ff(2, 5, 0.3) - 0.83692 + np.float64(0.83692) """ return binom.cdf(x, n, p) @@ -129,7 +129,7 @@ def sf(self, x, n, p): -------- >>> from surpyval import Binomial >>> Binomial.sf(2, 5, 0.3) - 0.16308 + np.float64(0.16308) """ return binom.sf(x, n, p) @@ -215,7 +215,7 @@ def qf(self, q, n, p): -------- >>> from surpyval import Binomial >>> Binomial.qf(0.5, 5, 0.3) - 1.0 + np.float64(1.0) """ return binom.ppf(q, n, p) @@ -289,7 +289,7 @@ def moment(self, m, n, p): -------- >>> from surpyval import Binomial >>> Binomial.moment(1, 5, 0.3) - 1.5 + np.float64(1.5) """ return binom.moment(m, n, p) @@ -302,7 +302,7 @@ def entropy(self, n, p): -------- >>> from surpyval import Binomial >>> Binomial.entropy(5, 0.3) - 1.413614855283445 + np.float64(1.413614855283445) """ return binom.entropy(n, p) @@ -414,7 +414,7 @@ def from_params(self, params): >>> from surpyval import Binomial >>> model = Binomial.from_params([5, 0.3]) >>> model.mean() - 1.5 + np.float64(1.5) """ params = np.atleast_1d(np.asarray(params, dtype=float)) diff --git a/surpyval/univariate/parametric/distributions/custom_distribution.py b/surpyval/univariate/parametric/distributions/custom_distribution.py index 745fbf9b..4a3daa03 100644 --- a/surpyval/univariate/parametric/distributions/custom_distribution.py +++ b/surpyval/univariate/parametric/distributions/custom_distribution.py @@ -43,16 +43,16 @@ class CustomDistribution(ParametricFitter): >>> name = 'Gompertz' >>> >>> def Hf(x, *params): - >>> return params[0] * np.exp(params[1] * x - 1) - >>> + ... return params[0] * np.exp(params[1] * x - 1) + ... >>> param_names = ['nu', 'b'] >>> bounds = ((0, None), (0, None)) >>> support = (-np.inf, np.inf) >>> Gompertz = surv.CustomDistribution( - name, Hf, param_names, bounds, support - ) + ... name, Hf, param_names, bounds, support + ... ) >>> x = np.array([1, 2, 3, 4, 5]) - >>> Gompertz.fit(x) + >>> model = Gompertz.fit(x) """ def __init__(self, name, fun, param_names, bounds, support): diff --git a/surpyval/univariate/parametric/distributions/expo_weibull.py b/surpyval/univariate/parametric/distributions/expo_weibull.py index 2cc80eee..bd7c337e 100755 --- a/surpyval/univariate/parametric/distributions/expo_weibull.py +++ b/surpyval/univariate/parametric/distributions/expo_weibull.py @@ -184,7 +184,7 @@ def cs(self, x, X, alpha, beta, mu): >>> import numpy as np >>> from surpyval import ExpoWeibull >>> x = np.array([1, 2, 3, 4, 5]) - >>> ExpoWeibull.sf(x, 1, 3, 4, 1.2) + >>> ExpoWeibull.cs(x, 1, 3, 4, 1.2) array([8.77367129e-01, 4.25451775e-01, 5.09266354e-02, 5.37452200e-04, 1.35732908e-07]) """ diff --git a/surpyval/univariate/parametric/distributions/exponential.py b/surpyval/univariate/parametric/distributions/exponential.py index 7ea07526..f760ebcd 100755 --- a/surpyval/univariate/parametric/distributions/exponential.py +++ b/surpyval/univariate/parametric/distributions/exponential.py @@ -327,7 +327,7 @@ def qf(self, p, failure_rate): >>> from surpyval import Exponential >>> p = np.array([.1, .2, .3, .4, .5]) >>> Exponential.qf(p, 3) - array([0.03512219, 0.07438118, 0.11889152, 0.17027853, 0.23104906]) + array([0.03512017, 0.07438118, 0.11889165, 0.17027521, 0.23104906]) """ return -np.log1p(-p) / failure_rate @@ -386,7 +386,7 @@ def moment(self, n, failure_rate): -------- >>> from surpyval import Exponential >>> Exponential.moment(2, 3) - 0.2222222222222222 + np.float64(0.2222222222222222) """ return factorial(n) / (failure_rate**n) @@ -414,7 +414,7 @@ def entropy(self, failure_rate): -------- >>> from surpyval import Exponential >>> Exponential.entropy(3) - -0.09861228866810978 + np.float64(-0.09861228866810978) """ return 1 - np.log(failure_rate) diff --git a/surpyval/univariate/parametric/distributions/gamma.py b/surpyval/univariate/parametric/distributions/gamma.py index 51a32c8e..2f7450d1 100755 --- a/surpyval/univariate/parametric/distributions/gamma.py +++ b/surpyval/univariate/parametric/distributions/gamma.py @@ -1,13 +1,8 @@ -import warnings - from autograd.scipy.special import gamma as agamma from autograd.scipy.special import gammaln as agammaln -from scipy.optimize import minimize from scipy.special import digamma, gammaincinv -from scipy.stats import pearsonr from surpyval import np -from surpyval.univariate.nonparametric import plotting_positions from surpyval.univariate.parametric.parametric_fitter import ParametricFitter from surpyval.utils.autograd_gamma_compat import gammainc as agammainc from surpyval.utils.autograd_gamma_compat import gammainccln as agammainccln @@ -35,6 +30,26 @@ def __init__(self, name): param_map={"alpha": 0, "beta": 1}, plot_x_scale="linear", ) + # The Gamma has no linearising probability plot, for the same + # reason as the Beta above it: the CDF is the regularised + # incomplete gamma function, and the shape sits *inside* that + # special function rather than outside it as an exponent. The + # only straight-line y-axis is the inverse incomplete gamma, + # which needs the shape -- so to draw the axis you need the + # answer, and to get the answer you need the axis. + # + # MPP broke the circle by guessing the shape from moments, + # drawing the plot on that guess and regressing. When the guess + # is off the axis is the wrong axis, the points are no longer + # straight on it, and the regression fits a line through a + # curve -- returning a confident, wrong estimate rather than an + # error. An offset makes it worse: the shift distorts the low-x + # end hardest, which is exactly where the shape information is. + # + # Fit by MLE (the default), MSE or MOM instead. ``plot()`` still + # works, because it transforms with the *fitted* parameters, so + # the axis is the right one by the time it is drawn. + self.supports_mpp = False @staticmethod def _moment_estimate(x): @@ -396,7 +411,7 @@ def moment(self, n, alpha, beta): -------- >>> from surpyval import Gamma >>> Gamma.moment(3, 3, 4) - 0.9375 + np.float64(0.9375) """ return agamma(n + alpha) / (beta**n * agamma(alpha)) @@ -430,7 +445,7 @@ def entropy(self, alpha, beta): -------- >>> from surpyval import Gamma >>> Gamma.entropy(3, 4) - 0.46128414924312033 + np.float64(0.46128414924312033) """ return ( alpha @@ -490,121 +505,5 @@ def mpp_inv_y_transform(self, y, *params): def mpp_x_transform(self, x, gamma=0): return x - gamma - def mpp( - self, - x, - c=None, - n=None, - t=None, - heuristic="Nelson-Aalen", - rr="y", - on_d_is_0=False, - offset=False, - ): - # Forward the truncation windows (previously dropped, #280). - x_pp, r, d, F = plotting_positions( - x, c=c, n=n, t=t, heuristic=heuristic - ) - - results = {} - - if on_d_is_0: - pass - else: - F = F[d > 0] - x_pp = x_pp[d > 0] - - if (F == 1).any(): - mask = F != 1 - warnings.warn( - "Some heuristic values for CDF = 1 have been " - "encountered in plotting points and have been " - "ignored.", - stacklevel=2, - ) - F = F[mask] - x_pp = x_pp[mask] - - init = self._parameter_initialiser(x_pp, c, n) - - mask = np.isfinite(F) - if not mask.all(): - warnings.warn( - "Some Infinite values encountered in plotting " - "points and have been ignored.", - stacklevel=2, - ) - F = F[mask] - x_pp = x_pp[mask] - - if offset: - - def fun(a): - return -pearsonr(x_pp, self.mpp_y_transform(F, a, 1.0))[0] - - # The moment-based init is computed from the *unshifted* data, - # which diverges for strongly offset data (tiny relative spread - # -> huge alpha) and strands the correlation search at a bad - # local optimum. Try several starts and keep the best (#257). - starts = {float(init[0]), 0.5, 1.0, 2.0, 5.0} - best = None - for a0 in starts: - res = minimize(fun, [a0], bounds=((1e-8, None),)) - if best is None or res.fun < best.fun: - best = res - alpha = best.x[0] - - y_pp = self.mpp_y_transform(F, alpha) - - if rr == "y": - # y = beta * (x - gamma): slope is beta, intercept is - # -beta * gamma. - params = np.polyfit(x_pp, y_pp, 1) - beta = params[0] - gamma = -params[1] / beta - elif rr == "x": - # x = y / beta + gamma: slope is 1/beta and the intercept - # IS gamma (#257). - params = np.polyfit(y_pp, x_pp, 1) - beta = 1.0 / params[0] - gamma = params[1] - - results["gamma"] = gamma - results["params"] = np.array([alpha, beta]) - - return results - else: - if rr == "y": - x_pp = x_pp[:, np.newaxis] - - def fun(alpha): - y_pp = self.mpp_y_transform(F, alpha, 1.0) - return np.linalg.lstsq(x_pp, y_pp)[1] - - res = minimize(fun, init[0], bounds=((0, None),)) - alpha = res.x[0] - y_pp = self.mpp_y_transform(F, alpha, 1.0) - beta, residuals, _, _ = np.linalg.lstsq(x_pp, y_pp) - beta = beta[0] - else: - # Regress against the same filtered plotting positions used - # everywhere else — the raw ``x`` argument has a different - # length whenever any point was filtered (#257). - - def fun(a): - y = self.mpp_y_transform(F, a, 1.0)[:, np.newaxis] - return np.linalg.lstsq(y, x_pp)[1] - - res = minimize(fun, init[0], bounds=((0, None),)) - alpha = res.x[0] - beta = np.linalg.lstsq( - self.mpp_y_transform(F, alpha, 1.0)[:, np.newaxis], x_pp - )[0][0] - beta = 1.0 / beta - - results["params"] = np.array([alpha, beta]) - - return results - Gamma: ParametricFitter = Gamma_("Gamma") diff --git a/surpyval/univariate/parametric/distributions/gumbel.py b/surpyval/univariate/parametric/distributions/gumbel.py index efd8aceb..4dc5bbd7 100755 --- a/surpyval/univariate/parametric/distributions/gumbel.py +++ b/surpyval/univariate/parametric/distributions/gumbel.py @@ -243,10 +243,9 @@ def qf(self, p, mu, sigma): -------- >>> import numpy as np >>> from surpyval import Gumbel - >>> p = np.array([.1, .2, .3, .4, .5]) + >>> p = np.array([0.1, 0.3, 0.5]) >>> Gumbel.qf(p, 3, 2) - array([-1.50073465e+00, 1.20026481e-04, 9.38139134e-01, 1.65654602e+00, - 2.26697416e+00]) + array([-1.50073465, 0.93813913, 2.26697416]) """ return mu + sigma * (np.log(-np.log1p(-p))) @@ -325,7 +324,7 @@ def entropy(self, mu, sigma): -------- >>> from surpyval import Gumbel >>> Gumbel.entropy(3, 2) - 2.270362845461478 + np.float64(2.270362845461478) """ return np.log(sigma) + euler_gamma + 1 diff --git a/surpyval/univariate/parametric/distributions/gumbel_lev.py b/surpyval/univariate/parametric/distributions/gumbel_lev.py index 352cfc76..8b5f3639 100644 --- a/surpyval/univariate/parametric/distributions/gumbel_lev.py +++ b/surpyval/univariate/parametric/distributions/gumbel_lev.py @@ -335,7 +335,7 @@ def entropy(self, mu, sigma): -------- >>> from surpyval import GumbelLEV >>> GumbelLEV.entropy(3, 2) - 2.270362845461478 + np.float64(2.270362845461478) """ return np.log(sigma) + euler_gamma + 1 diff --git a/surpyval/univariate/parametric/distributions/logistic.py b/surpyval/univariate/parametric/distributions/logistic.py index 88924a63..53cdd424 100755 --- a/surpyval/univariate/parametric/distributions/logistic.py +++ b/surpyval/univariate/parametric/distributions/logistic.py @@ -226,9 +226,9 @@ def qf(self, p, mu, sigma): -------- >>> import numpy as np >>> from surpyval import Logistic - >>> p = np.array([.1, .2, .3, .4, .5]) + >>> p = np.array([0.1, 0.2, 0.3, 0.4]) >>> Logistic.qf(p, 3, 4) - array([-5.78889831, -2.54517744, -0.38919144, 1.37813957, 3. ]) + array([-5.78889831, -2.54517744, -0.38919144, 1.37813957]) """ return mu + sigma * (np.log(p) - np.log1p(-p)) @@ -311,7 +311,7 @@ def entropy(self, mu, sigma): -------- >>> from surpyval import Logistic >>> Logistic.entropy(3, 4) - 3.386294361119891 + np.float64(3.386294361119891) """ return np.log(sigma) + 2 diff --git a/surpyval/univariate/parametric/distributions/loglogistic.py b/surpyval/univariate/parametric/distributions/loglogistic.py index 597b6dca..54b0b215 100755 --- a/surpyval/univariate/parametric/distributions/loglogistic.py +++ b/surpyval/univariate/parametric/distributions/loglogistic.py @@ -57,7 +57,7 @@ def sf(self, x, alpha, beta): >>> from surpyval import LogLogistic >>> x = np.array([1, 2, 3, 4, 5]) >>> LogLogistic.sf(x, 3, 4) - array([0.62245933, 0.5621765 , 0.5 , 0.4378235 , 0.37754067]) + array([0.98780488, 0.83505155, 0.5 , 0.24035608, 0.11473088]) """ # 1 / (1 + (x/alpha)^beta): algebraically identical to the # (x/alpha)^-beta form but defined at x = 0 (sf(0) = 1) instead @@ -304,7 +304,7 @@ def mean(self, alpha, beta): -------- >>> from surpyval import LogLogistic >>> LogLogistic.mean(3, 4) - 3 + np.float64(3.332162203618775) """ if beta > 1: return (alpha * np.pi / beta) / (np.sin(np.pi / beta)) @@ -382,7 +382,7 @@ def entropy(self, alpha, beta): -------- >>> from surpyval import LogLogistic >>> LogLogistic.entropy(3, 4) - 1.7123179275482192 + np.float64(1.7123179275482192) """ return np.log(alpha / beta) + 2 diff --git a/surpyval/univariate/parametric/distributions/lognormal.py b/surpyval/univariate/parametric/distributions/lognormal.py index 36ff05d4..571059a2 100755 --- a/surpyval/univariate/parametric/distributions/lognormal.py +++ b/surpyval/univariate/parametric/distributions/lognormal.py @@ -301,10 +301,9 @@ def qf(self, p, mu, sigma): -------- >>> import numpy as np >>> from surpyval import LogNormal - >>> p = np.array([.1, .2, .3, .4, .5]) + >>> p = np.array([0.1, 0.2, 0.3, 0.4]) >>> LogNormal.qf(p, 3, 4) - array([ 0.11928899, 0.69316658, 2.46550819, 7.29078766, - 20.08553692]) + array([0.11928899, 0.69316658, 2.46550819, 7.29078766]) """ return np.exp(scipy_norm.ppf(p, mu, sigma)) @@ -334,7 +333,7 @@ def mean(self, mu, sigma): -------- >>> from surpyval import LogNormal >>> LogNormal.mean(3, 4) - 59874.14171519782 + np.float64(59874.14171519782) """ return np.exp(mu + (sigma**2) / 2) @@ -366,7 +365,7 @@ def moment(self, n, mu, sigma): -------- >>> from surpyval import LogNormal >>> LogNormal.moment(2, 3, 4) - 3.1855931757113756e+16 + np.float64(3.1855931757113756e+16) """ return np.exp(n * mu + (n**2 * sigma**2) / 2) @@ -396,7 +395,7 @@ def entropy(self, mu, sigma): -------- >>> from surpyval import LogNormal >>> LogNormal.entropy(3, 4) - 5.805232894324563 + np.float64(5.805232894324563) """ return mu + 0.5 * np.log(2 * np.pi * np.e * sigma**2) diff --git a/surpyval/univariate/parametric/distributions/normal.py b/surpyval/univariate/parametric/distributions/normal.py index 3afcd536..8b34429c 100755 --- a/surpyval/univariate/parametric/distributions/normal.py +++ b/surpyval/univariate/parametric/distributions/normal.py @@ -308,9 +308,9 @@ def qf(self, p, mu, sigma): -------- >>> import numpy as np >>> from surpyval import Normal - >>> p = np.array([.1, .2, .3, .4, .5]) + >>> p = np.array([0.1, 0.2, 0.3, 0.4]) >>> Normal.qf(p, 3, 4) - array([-2.12620626, -0.36648493, 0.90239795, 1.98661159, 3. ]) + array([-2.12620626, -0.36648493, 0.90239795, 1.98661159]) """ return scipy_norm.ppf(p, mu, sigma) @@ -372,7 +372,7 @@ def moment(self, n, mu, sigma): -------- >>> from surpyval import Normal >>> Normal.moment(2, 3, 4) - 25.0 + np.float64(25.0) """ return scipy_norm.moment(n, mu, sigma) @@ -402,7 +402,7 @@ def entropy(self, mu, sigma): -------- >>> from surpyval import Normal >>> Normal.entropy(3, 4) - 2.8052328943245635 + np.float64(2.8052328943245635) """ return 0.5 * np.log(2 * np.pi * np.e * sigma**2) diff --git a/surpyval/univariate/parametric/distributions/rayleigh.py b/surpyval/univariate/parametric/distributions/rayleigh.py index 98e72f54..719ac981 100644 --- a/surpyval/univariate/parametric/distributions/rayleigh.py +++ b/surpyval/univariate/parametric/distributions/rayleigh.py @@ -299,7 +299,7 @@ def mean(self, sigma): -------- >>> from surpyval import Rayleigh >>> Rayleigh.mean(3) - 3.7599424119465006 + np.float64(3.7599424119465006) """ return sigma * np.sqrt(np.pi / 2) @@ -329,7 +329,7 @@ def moment(self, n, sigma): -------- >>> from surpyval import Rayleigh >>> Rayleigh.moment(2, 3) - 18.0 + np.float64(18.0) """ return (sigma**n) * (2 ** (n / 2)) * gamma_func(1 + n / 2) diff --git a/surpyval/univariate/parametric/distributions/uniform.py b/surpyval/univariate/parametric/distributions/uniform.py index e5fb9dd0..608eb615 100755 --- a/surpyval/univariate/parametric/distributions/uniform.py +++ b/surpyval/univariate/parametric/distributions/uniform.py @@ -119,7 +119,7 @@ def ff(self, x, a, b): >>> import numpy as np >>> from surpyval import Uniform >>> x = np.array([1, 2, 3, 4, 5]) - >>> Uniform.sf(x, 0, 6) + >>> Uniform.ff(x, 0, 6) array([0.16666667, 0.33333333, 0.5 , 0.66666667, 0.83333333]) """ f = np.zeros_like(x) @@ -341,7 +341,7 @@ def moment(self, n, a, b): -------- >>> from surpyval import Uniform >>> Uniform.moment(2, 0, 6) - 12.0 + np.float64(12.0) """ if n == 0: return 1 @@ -377,7 +377,7 @@ def entropy(self, a, b): -------- >>> from surpyval import Uniform >>> Uniform.entropy(0, 6) - 1.791759469228055 + np.float64(1.791759469228055) """ return np.log(b - a) diff --git a/surpyval/univariate/parametric/distributions/weibull.py b/surpyval/univariate/parametric/distributions/weibull.py index 59e8ac5a..f66a6bd0 100755 --- a/surpyval/univariate/parametric/distributions/weibull.py +++ b/surpyval/univariate/parametric/distributions/weibull.py @@ -302,7 +302,7 @@ def mean(self, alpha, beta): -------- >>> from surpyval import Weibull >>> Weibull.mean(3, 4) - 2.7192074311664314 + np.float64(2.7192074311664314) """ return alpha * gamma_func(1 + 1.0 / beta) @@ -334,7 +334,7 @@ def moment(self, n, alpha, beta): -------- >>> from surpyval import Weibull >>> Weibull.moment(2, 3, 4) - 7.976042329074821 + np.float64(7.976042329074821) """ return alpha**n * gamma_func(1 + n / beta) diff --git a/surpyval/univariate/parametric/mixture_model.py b/surpyval/univariate/parametric/mixture_model.py index 6615d657..2a6f80ea 100755 --- a/surpyval/univariate/parametric/mixture_model.py +++ b/surpyval/univariate/parametric/mixture_model.py @@ -297,11 +297,11 @@ def fit( Sub-Distributions : 2 Fitted by : EM Weights : - 0.6094710980384728, - 0.39052890196152723 + 0.6184891886499861, + 0.381510811350014 Parameters : - alpha: [ 5.8855232 17.23187124] - beta: [ 2.04051304 11.01565277] + alpha: [ 6.32508961 17.37701969] + beta: [ 1.83105154 12.01392721] """ data = SurpyvalData(x=x, c=c, n=n, t=t, tl=tl, tr=tr, xl=xl, xr=xr) diff --git a/surpyval/univariate/parametric/parametric.py b/surpyval/univariate/parametric/parametric.py index 9e7a25c6..a0b1dcb7 100755 --- a/surpyval/univariate/parametric/parametric.py +++ b/surpyval/univariate/parametric/parametric.py @@ -549,7 +549,7 @@ def sf(self, x: npt.ArrayLike) -> npt.NDArray: >>> from surpyval import Weibull >>> model = Weibull.from_params([10, 3]) >>> model.sf(2) - 0.9920319148370607 + np.float64(0.9920319148370607) >>> model.sf([1, 2, 3, 4, 5]) array([0.9990005 , 0.99203191, 0.97336124, 0.938005 , 0.8824969 ]) """ @@ -591,7 +591,7 @@ def ff(self, x: npt.ArrayLike) -> npt.NDArray: >>> from surpyval import Weibull >>> model = Weibull.from_params([10, 3]) >>> model.ff(2) - 0.007968085162939342 + np.float64(0.007968085162939372) >>> model.ff([1, 2, 3, 4, 5]) array([0.0009995 , 0.00796809, 0.02663876, 0.061995 , 0.1175031 ]) """ @@ -632,7 +632,7 @@ def df(self, x: npt.ArrayLike) -> npt.NDArray: >>> from surpyval import Weibull >>> model = Weibull.from_params([10, 3]) >>> model.df(2) - 0.01190438297804473 + np.float64(0.01190438297804473) >>> model.df([1, 2, 3, 4, 5]) array([0.002997 , 0.01190438, 0.02628075, 0.04502424, 0.06618727]) """ @@ -679,7 +679,7 @@ def hf(self, x: npt.ArrayLike) -> npt.NDArray: >>> from surpyval import Weibull >>> model = Weibull.from_params([10, 3]) >>> model.hf(2) - 0.012000000000000002 + np.float64(0.012000000000000002) >>> model.hf([1, 2, 3, 4, 5]) array([0.003, 0.012, 0.027, 0.048, 0.075]) """ @@ -687,7 +687,12 @@ def hf(self, x: npt.ArrayLike) -> npt.NDArray: if (self.p == 1) and (self.f0 == 0): xg = x - self.gamma # type: ignore[operator] s0 = getattr(self.dist, "support", (-np.inf, np.inf))[0] - return np.where(xg < s0, 0.0, self.dist.hf(xg, *self.params)) + out = np.where(xg < s0, 0.0, self.dist.hf(xg, *self.params)) + # ``np.where`` hands back a 0-d array for scalar input, so a + # scalar argument used to come out as ``array(0.012)`` while + # every sibling method returned a numpy scalar. ``[()]`` is a + # no-op on a real array and unwraps the 0-d case. + return out[()] else: return self.df(x) / self.sf(x) @@ -719,7 +724,7 @@ def Hf(self, x: npt.ArrayLike) -> npt.NDArray: >>> from surpyval import Weibull >>> model = Weibull.from_params([10, 3]) >>> model.Hf(2) - 0.008000000000000002 + np.float64(0.008000000000000002) >>> model.Hf([1, 2, 3, 4, 5]) array([0.001, 0.008, 0.027, 0.064, 0.125]) """ @@ -728,7 +733,8 @@ def Hf(self, x: npt.ArrayLike) -> npt.NDArray: if (self.p == 1) and (self.f0 == 0): xg = x - self.gamma # type: ignore[operator] s0 = getattr(self.dist, "support", (-np.inf, np.inf))[0] - return np.where(xg < s0, 0.0, self.dist.Hf(xg, *self.params)) + out = np.where(xg < s0, 0.0, self.dist.Hf(xg, *self.params)) + return out[()] else: return -np.log(self.sf(x)) @@ -757,7 +763,7 @@ def qf(self, p: npt.ArrayLike) -> npt.NDArray: >>> from surpyval import Weibull >>> model = Weibull.from_params([10, 3]) >>> model.qf(0.2) - 6.06542793124108 + np.float64(6.06542793124108) >>> model.qf([.1, .2, .3, .4, .5]) array([4.72308719, 6.06542793, 7.09181722, 7.99387877, 8.84997045]) @@ -823,14 +829,14 @@ def cs(self, x: npt.ArrayLike, X: npt.ArrayLike) -> npt.NDArray: >>> from surpyval import Weibull >>> model = Weibull.from_params([10, 3]) >>> model.cs(11, 10) - 0.00025840046151723767 + np.float64(0.00025840046151723767) """ x = np.asarray(x) X = np.asarray(X) Xg = X - self.gamma # type: ignore[operator] cs = np.array(self.dist.cs(x, Xg, *self.params)) cs[cs > 1.0] = 1 - return cs + return cs[()] def random( self, @@ -972,7 +978,7 @@ def mean(self) -> float: >>> from surpyval import Weibull >>> model = Weibull.from_params([10, 3]) >>> model.mean() - 8.929795115692489 + np.float64(8.929795115692489) """ if not hasattr(self, "_mean"): # Defective mean: the zero-inflated mass f0 sits at 0 and @@ -998,7 +1004,7 @@ def var(self) -> float: >>> from surpyval import Weibull >>> model = Weibull.from_params([10, 3]) >>> model.var() - 11.229... + np.float64(10.533288486847923) """ m1 = self.dist._moment(1, *self.params) m2 = self.dist._moment(2, *self.params) @@ -1071,7 +1077,7 @@ def entropy(self) -> float: >>> from surpyval import Normal >>> model = Normal.from_params([10, 3]) >>> model.entropy() - 2.5175508218727822 + np.float64(2.5175508218727822) Notes ----- @@ -1556,16 +1562,19 @@ def plot( Returns ------- - plot : list - list of a matplotlib plot object + plot : matplotlib.axes.Axes + the axes the probability plot was drawn onto Examples -------- + >>> import numpy as np >>> from surpyval import Weibull + >>> np.random.seed(1) >>> x = Weibull.random(100, 10, 3) >>> model = Weibull.fit(x) >>> model.plot() + """ if ax is None: ax = plt.gcf().gca() diff --git a/surpyval/univariate/parametric/parametric_fitter.py b/surpyval/univariate/parametric/parametric_fitter.py index 43d2842f..3c52b9bc 100755 --- a/surpyval/univariate/parametric/parametric_fitter.py +++ b/surpyval/univariate/parametric/parametric_fitter.py @@ -753,6 +753,7 @@ def fit( -------- >>> from surpyval import Weibull >>> import numpy as np + >>> np.random.seed(1) >>> x = Weibull.random(100, 10, 4) >>> model = Weibull.fit(x) >>> print(model) @@ -761,8 +762,8 @@ def fit( Distribution : Weibull Fitted by : MLE Parameters : - alpha: 10.551521182640098 - beta: 3.792549834495306 + alpha: 9.815018791049368 + beta: 3.798740470368033 >>> Weibull.fit(x, how='MPS', fixed={'alpha' : 10}) Parametric SurPyval Model ========================= @@ -770,15 +771,16 @@ def fit( Fitted by : MPS Parameters : alpha: 10.0 - beta: 3.4314657446866836 - >>> Weibull.fit(xl=x-1, xr=x+1, how='MPP') + beta: 3.670796510564323 + >>> Weibull.fit(xl=np.floor(x), xr=np.ceil(x), how='MPP', + ... heuristic='Turnbull') Parametric SurPyval Model ========================= Distribution : Weibull Fitted by : MPP Parameters : - alpha: 9.943092756713078 - beta: 8.613016934518258 + alpha: 9.950168329892755 + beta: 3.211971411540382 >>> c = np.zeros_like(x) >>> c[x > 13] = 1 >>> x[x > 13] = 13 @@ -790,8 +792,8 @@ def fit( Distribution : Weibull Fitted by : MLE Parameters : - alpha: 10.363725328793413 - beta: 4.9886821457305865 + alpha: 9.893584496413128 + beta: 3.78688602908912 """ surv_data = SurpyvalData( @@ -892,10 +894,10 @@ def fit_from_df( ========================= Distribution : Weibull Fitted by : MLE - Offset (gamma) : 39.76562962867477 + Offset (gamma) : 39.76557772434183 Parameters : - alpha: 7.141925216146524 - beta: 2.6204524040137844 + alpha: 7.141983615103902 + beta: 2.62047590823775 """ if not isinstance(df, pd.DataFrame): diff --git a/surpyval/univariate/regression/accelerated_failure_time/aft_fitter.py b/surpyval/univariate/regression/accelerated_failure_time/aft_fitter.py index 13d594e3..ae87b354 100644 --- a/surpyval/univariate/regression/accelerated_failure_time/aft_fitter.py +++ b/surpyval/univariate/regression/accelerated_failure_time/aft_fitter.py @@ -144,8 +144,14 @@ def AFT(distribution): Examples -------- + >>> import numpy as np >>> from surpyval import Weibull >>> from surpyval import AFT - >>> model = AFT(Weibull).fit(x, Z=covariates, c=c) + >>> np.random.seed(1) + >>> Z = np.random.binomial(1, 0.5, 100).reshape(-1, 1) + >>> x = Weibull.random(100, 10, 2) * np.exp(-0.5 * Z[:, 0]) + >>> model = AFT(Weibull).fit(x, Z=Z) + >>> model.params.round(3) + array([9.629, 1.751, 0.473]) """ return AFTFitter(distribution) diff --git a/surpyval/univariate/regression/accelerated_life/accelerated_life.py b/surpyval/univariate/regression/accelerated_life/accelerated_life.py index 3fa5fc39..4626f273 100644 --- a/surpyval/univariate/regression/accelerated_life/accelerated_life.py +++ b/surpyval/univariate/regression/accelerated_life/accelerated_life.py @@ -38,9 +38,15 @@ def AcceleratedLife(distribution, life_model): Examples -------- + >>> import numpy as np >>> from surpyval import Weibull >>> from surpyval import AcceleratedLife, Power - >>> model = AcceleratedLife(Weibull, Power).fit(x, c=c, Z=stress) + >>> np.random.seed(1) + >>> stress = np.repeat([20.0, 30.0, 40.0], 40).reshape(-1, 1) + >>> x = Weibull.random(120, 10, 3) * (100.0 / stress[:, 0]) + >>> model = AcceleratedLife(Weibull, Power).fit(x, Z=stress) + >>> model.params.round(3) + array([ 1. , 2.831, 558.686, -0.828]) """ if distribution.name not in _LIFE_PARAM_MAP: supported = list(_LIFE_PARAM_MAP.keys()) diff --git a/surpyval/univariate/regression/additive_hazards/__init__.py b/surpyval/univariate/regression/additive_hazards/__init__.py index da90dd75..9a840da1 100644 --- a/surpyval/univariate/regression/additive_hazards/__init__.py +++ b/surpyval/univariate/regression/additive_hazards/__init__.py @@ -32,8 +32,14 @@ def AH(distribution): Examples -------- + >>> import numpy as np >>> from surpyval import Weibull, AH - >>> model = AH(Weibull).fit(x, Z=covariates, c=c) + >>> np.random.seed(1) + >>> Z = np.random.binomial(1, 0.5, 100).reshape(-1, 1) + >>> x = Weibull.random(100, 10, 2) * np.exp(-0.5 * Z[:, 0]) + >>> model = AH(Weibull).fit(x, Z=Z) + >>> model.params.round(3) + array([9.332, 1.851, 0.086]) """ return AdditiveHazardsFitter.create(distribution) diff --git a/surpyval/univariate/regression/additive_hazards/additive_hazards_fitter.py b/surpyval/univariate/regression/additive_hazards/additive_hazards_fitter.py index 59c7c13e..b2115f4a 100644 --- a/surpyval/univariate/regression/additive_hazards/additive_hazards_fitter.py +++ b/surpyval/univariate/regression/additive_hazards/additive_hazards_fitter.py @@ -200,8 +200,15 @@ def fit( Examples -------- - >>> from surpyval import WeibullAH + >>> import numpy as np + >>> from surpyval import Weibull, WeibullAH + >>> np.random.seed(1) + >>> Z = np.random.binomial(1, 0.5, 100).reshape(-1, 1) + >>> x = Weibull.random(100, 10, 2) * np.exp(-0.5 * Z[:, 0]) + >>> c = np.zeros(100) >>> model = WeibullAH.fit(x=x, Z=Z, c=c) + >>> model.params.round(3) + array([9.332, 1.851, 0.086]) """ data, prep = prepare_regression_fit( self, diff --git a/surpyval/univariate/regression/frailty/__init__.py b/surpyval/univariate/regression/frailty/__init__.py index c09dfeaa..724df8bd 100644 --- a/surpyval/univariate/regression/frailty/__init__.py +++ b/surpyval/univariate/regression/frailty/__init__.py @@ -34,8 +34,17 @@ def Frailty(distribution, family="gamma"): Examples -------- + >>> import numpy as np >>> from surpyval import Frailty, Weibull - >>> model = Frailty(Weibull).fit(x, Z=Z, c=c, groups=unit_id) + >>> np.random.seed(1) + >>> unit_id = np.repeat(np.arange(20), 5) + >>> Z = np.random.binomial(1, 0.5, 100).reshape(-1, 1) + >>> x = Weibull.random(100, 10, 2) * np.exp(-0.5 * Z[:, 0]) + >>> model = Frailty(Weibull).fit(x, Z=Z, groups=unit_id) + >>> model.beta.round(3) + array([0.829]) + >>> model.n_groups + 20 """ return FrailtyFitter.create(distribution, family) diff --git a/surpyval/univariate/regression/parametric_regression_model.py b/surpyval/univariate/regression/parametric_regression_model.py index a36c8664..efa9d8e7 100644 --- a/surpyval/univariate/regression/parametric_regression_model.py +++ b/surpyval/univariate/regression/parametric_regression_model.py @@ -430,12 +430,14 @@ def sf( Examples -------- - >>> from surpyval import Weibull - >>> model = Weibull.from_params([10, 3]) - >>> model.sf(2) - 0.9920319148370607 - >>> model.sf([1, 2, 3, 4, 5]) - array([0.9990005 , 0.99203191, 0.97336124, 0.938005 , 0.8824969 ]) + >>> import numpy as np + >>> from surpyval import Weibull, WeibullPH + >>> np.random.seed(1) + >>> Z = np.random.binomial(1, 0.5, 100).reshape(-1, 1) + >>> x = Weibull.random(100, 10, 2) * np.exp(-0.5 * Z[:, 0]) + >>> model = WeibullPH.fit(x, Z) + >>> model.sf([1, 2, 3], [[0], [0], [1]]).round(4) + array([0.9812, 0.9382, 0.7429]) """ if isinstance(x, list): x = np.array(x) @@ -669,12 +671,14 @@ def ff( Examples -------- - >>> from surpyval import Weibull - >>> model = Weibull.from_params([10, 3]) - >>> model.ff(2) - 0.007968085162939342 - >>> model.ff([1, 2, 3, 4, 5]) - array([0.0009995 , 0.00796809, 0.02663876, 0.061995 , 0.1175031 ]) + >>> import numpy as np + >>> from surpyval import Weibull, WeibullPH + >>> np.random.seed(1) + >>> Z = np.random.binomial(1, 0.5, 100).reshape(-1, 1) + >>> x = Weibull.random(100, 10, 2) * np.exp(-0.5 * Z[:, 0]) + >>> model = WeibullPH.fit(x, Z) + >>> model.ff([1, 2, 3], [[0], [0], [1]]).round(4) + array([0.0188, 0.0618, 0.2571]) """ if isinstance(x, list): x = np.array(x) @@ -708,12 +712,14 @@ def df( Examples -------- - >>> from surpyval import Weibull - >>> model = Weibull.from_params([10, 3]) - >>> model.df(2) - 0.01190438297804473 - >>> model.df([1, 2, 3, 4, 5]) - array([0.002997 , 0.01190438, 0.02628075, 0.04502424, 0.06618727]) + >>> import numpy as np + >>> from surpyval import Weibull, WeibullPH + >>> np.random.seed(1) + >>> Z = np.random.binomial(1, 0.5, 100).reshape(-1, 1) + >>> x = Weibull.random(100, 10, 2) * np.exp(-0.5 * Z[:, 0]) + >>> model = WeibullPH.fit(x, Z) + >>> model.df([1, 2, 3], [[0], [0], [1]]).round(4) + array([0.0326, 0.0524, 0.1289]) """ if isinstance(x, list): x = np.array(x) @@ -748,12 +754,14 @@ def hf( Examples -------- - >>> from surpyval import Weibull - >>> model = Weibull.from_params([10, 3]) - >>> model.hf(2) - 0.012000000000000002 - >>> model.hf([1, 2, 3, 4, 5]) - array([0.003, 0.012, 0.027, 0.048, 0.075]) + >>> import numpy as np + >>> from surpyval import Weibull, WeibullPH + >>> np.random.seed(1) + >>> Z = np.random.binomial(1, 0.5, 100).reshape(-1, 1) + >>> x = Weibull.random(100, 10, 2) * np.exp(-0.5 * Z[:, 0]) + >>> model = WeibullPH.fit(x, Z) + >>> model.hf([1, 2, 3], [[0], [0], [1]]).round(4) + array([0.0332, 0.0559, 0.1735]) """ if isinstance(x, list): x = np.array(x) @@ -789,12 +797,14 @@ def Hf( Examples -------- - >>> from surpyval import Weibull - >>> model = Weibull.from_params([10, 3]) - >>> model.Hf(2) - 0.008000000000000002 - >>> model.Hf([1, 2, 3, 4, 5]) - array([0.001, 0.008, 0.027, 0.064, 0.125]) + >>> import numpy as np + >>> from surpyval import Weibull, WeibullPH + >>> np.random.seed(1) + >>> Z = np.random.binomial(1, 0.5, 100).reshape(-1, 1) + >>> x = Weibull.random(100, 10, 2) * np.exp(-0.5 * Z[:, 0]) + >>> model = WeibullPH.fit(x, Z) + >>> model.Hf([1, 2, 3], [[0], [0], [1]]).round(4) + array([0.0189, 0.0638, 0.2972]) """ if isinstance(x, list): x = np.array(x) @@ -826,15 +836,22 @@ def random( Examples -------- - >>> from surpyval import Weibull - >>> model = Weibull.from_params([10, 3]) - >>> np.random.seed(1) - >>> model.random(1) - array([8.14127103]) - >>> from surpyval import WeibullPH + >>> import numpy as np + >>> from surpyval import Weibull, WeibullPH >>> np.random.seed(1) + >>> Z = np.random.binomial(1, 0.5, 100).reshape(-1, 1) + >>> x = Weibull.random(100, 10, 2) * np.exp(-0.5 * Z[:, 0]) >>> model = WeibullPH.fit(x, Z) - >>> x_rand, Z_rand = model.random(10, Z[:1]) + >>> np.random.seed(1) + >>> x_rand, Z_rand = model.random(5, Z[:1]) + >>> x_rand.round(3) + array([ 8.919, 5.095, 33.929, 10.666, 13.97 ]) + >>> Z_rand + array([[0.], + [0.], + [0.], + [0.], + [0.]]) """ # Dispatch to the regression fitter's own covariate-aware sampler # (#261): the previous implementation ignored ``Z`` entirely and diff --git a/surpyval/univariate/regression/proportional_hazards/__init__.py b/surpyval/univariate/regression/proportional_hazards/__init__.py index ae74e0fa..add16b8c 100644 --- a/surpyval/univariate/regression/proportional_hazards/__init__.py +++ b/surpyval/univariate/regression/proportional_hazards/__init__.py @@ -31,9 +31,15 @@ def PH(distribution): Examples -------- + >>> import numpy as np >>> from surpyval import Weibull >>> from surpyval import PH - >>> model = PH(Weibull).fit(x, Z=covariates, c=c) + >>> np.random.seed(1) + >>> Z = np.random.binomial(1, 0.5, 100).reshape(-1, 1) + >>> x = Weibull.random(100, 10, 2) * np.exp(-0.5 * Z[:, 0]) + >>> model = PH(Weibull).fit(x, Z=Z) + >>> model.params.round(3) + array([9.629, 1.751, 0.829]) """ return ProportionalHazardsFitter.create(distribution) diff --git a/surpyval/univariate/regression/proportional_hazards/proportional_hazards_fitter.py b/surpyval/univariate/regression/proportional_hazards/proportional_hazards_fitter.py index f0f17f64..8497d0bc 100644 --- a/surpyval/univariate/regression/proportional_hazards/proportional_hazards_fitter.py +++ b/surpyval/univariate/regression/proportional_hazards/proportional_hazards_fitter.py @@ -188,17 +188,13 @@ def fit( >>> from surpyval import WeibullPH >>> from surpyval.datasets import load_tires_data - >>> from autograd import numpy as anp - >>> import numpy as np - >>> >>> data = load_tires_data() - >>> >>> x = data['Survival'].values >>> c = data['Censoring'].values >>> Z = data[[ - 'Wedge gauge', 'Interbelt gauge', 'Peel force', - 'Wedge gauge×peel force' - ]].values + ... 'Wedge gauge', 'Interbelt gauge', 'Peel force', + ... 'Wedge gauge×peel force' + ... ]].values >>> model = WeibullPH.fit(x=x, Z=Z, c=c) >>> model Parametric Regression SurPyval Model @@ -208,13 +204,13 @@ def fit( Regression Model : Log Linear [e^(beta'Z)] Fitted by : MLE Distribution : - alpha: 0.24255054642143947 - beta: 16.057791674515805 + alpha: 0.2425513627560218 + beta: 16.057785182711932 Regression Model : - beta_0: -9.165062641226692 - beta_1: -7.998599877425742 - beta_2: -27.503283340963034 - beta_3: 18.38550143851751 + beta_0: -9.165062726518311 + beta_1: -7.998573055929788 + beta_2: -27.50318580568538 + beta_3: 18.385445332039488 >>> model = WeibullPH.fit(x=x, Z=Z, c=c, fixed={"beta": 15}) >>> model Parametric Regression SurPyval Model @@ -224,13 +220,13 @@ def fit( Regression Model : Log Linear [e^(beta'Z)] Fitted by : MLE Distribution : - alpha: 0.23772915681951018 - beta: 15.0 + alpha: 0.237729668424067 + beta: 15.0 Regression Model : - beta_0: -8.628333861229965 - beta_1: -7.617541980158942 - beta_2: -25.952407717383302 - beta_3: 17.270173771235655 + beta_0: -8.62832691738283 + beta_1: -7.617529362323243 + beta_2: -25.952367249502934 + beta_3: 17.270148387391387 """ data, prep = prepare_regression_fit( self, diff --git a/surpyval/univariate/regression/proportional_odds/proportional_odds_fitter.py b/surpyval/univariate/regression/proportional_odds/proportional_odds_fitter.py index e6074230..d0304bbd 100644 --- a/surpyval/univariate/regression/proportional_odds/proportional_odds_fitter.py +++ b/surpyval/univariate/regression/proportional_odds/proportional_odds_fitter.py @@ -200,8 +200,14 @@ def PO(distribution): Examples -------- + >>> import numpy as np >>> from surpyval import Logistic >>> from surpyval import PO - >>> model = PO(Logistic).fit(x, Z=covariates, c=c) + >>> np.random.seed(1) + >>> Z = np.random.binomial(1, 0.5, 100).reshape(-1, 1) + >>> x = Logistic.random(100, 10, 2) + 2.0 * Z[:, 0] + >>> model = PO(Logistic).fit(x, Z=Z) + >>> model.params.round(3) + array([9.708, 2.337, 0.918]) """ return ProportionalOddsFitter(distribution) diff --git a/surpyval/univariate/regression/regression_data.py b/surpyval/univariate/regression/regression_data.py index 259390ec..b85bf809 100644 --- a/surpyval/univariate/regression/regression_data.py +++ b/surpyval/univariate/regression/regression_data.py @@ -307,11 +307,26 @@ def fit_from_df( Examples -------- - >>> from surpyval import WeibullPH + >>> import numpy as np + >>> import pandas as pd + >>> from surpyval import Weibull, WeibullPH + >>> np.random.seed(1) + >>> age = np.random.uniform(20, 60, 100) + >>> weight = np.random.uniform(50, 100, 100) + >>> time = Weibull.random(100, 10, 2) * np.exp(-0.02 * (age - 40)) + >>> df = pd.DataFrame({ + ... "time": time, + ... "age": age, + ... "weight": weight, + ... "censored": np.zeros(100, dtype=int), + ... }) >>> model = WeibullPH.fit_from_df( ... df, x_col="time", Z_cols=["age", "weight"], c_col="censored" ... ) - >>> model.sf([10, 20], df[["age", "weight"]]) + >>> model.feature_names + ['age', 'weight'] + >>> model.sf([10, 20], df[["age", "weight"]].head(2)).round(4) + array([0.4757, 0.0024]) """ Z, feature_names, model_spec = design_matrix_from_df( df, Z_cols, formula diff --git a/surpyval/univariate/regression/tvc_schedule.py b/surpyval/univariate/regression/tvc_schedule.py index f2589a4c..25ae5d8b 100644 --- a/surpyval/univariate/regression/tvc_schedule.py +++ b/surpyval/univariate/regression/tvc_schedule.py @@ -448,7 +448,9 @@ def from_expression(cls, expr, horizon, resolution=1.0, t0=0.0): Examples -------- >>> StepSchedule.from_expression("0.9 if t % 24 < 8 else 0.3", 96) + StepSchedule(step, 9 segment(s), p=1) >>> StepSchedule.from_expression("0.3 * 2 ** floor(t / 1000)", 5000) + StepSchedule(step, 6 segment(s), p=1) """ exprs = [expr] if isinstance(expr, str) else list(expr) if len(exprs) == 0: diff --git a/surpyval/utils/__init__.py b/surpyval/utils/__init__.py index 5f6064f9..604ae09a 100755 --- a/surpyval/utils/__init__.py +++ b/surpyval/utils/__init__.py @@ -302,7 +302,7 @@ def xrd_handler(x, r, d): >>> r array([5, 4, 3, 2, 1]) >>> d - array([1, 1, 1, 1, 1])) + array([1, 1, 1, 1, 1]) """ try: @@ -824,7 +824,7 @@ def xcnt_to_xrd(x, c=None, n=None, t=None, **kwargs): >>> n = np.array([1, 1, 1, 1, 1]) >>> x, r, d = xcnt_to_xrd(x, c, n) >>> x - array([1, 2, 3, 4, 5]) + array([1., 2., 3., 4., 5.]) >>> r array([5, 4, 3, 2, 1]) >>> d @@ -917,8 +917,8 @@ def xrd_to_xcnt(x, r, d): >>> r = np.array([5, 4, 3, 2, 1]) >>> d = np.array([1, 0, 0, 1, 1]) >>> x, c, n, t = xrd_to_xcnt(x, r, d) - >>> x, c, n, t - array([1, 2, 3, 4, 5]) + >>> x + array([1., 2., 3., 4., 5.]) >>> c array([0, 1, 1, 0, 0]) >>> n @@ -928,7 +928,7 @@ def xrd_to_xcnt(x, r, d): [-inf, inf], [-inf, inf], [-inf, inf], - [-inf, inf]])) + [-inf, inf]]) """ n_f = np.copy(d) x_f = np.copy(x) @@ -1001,7 +1001,7 @@ def fsli_to_xcnt(f=None, s=None, l=None, i=None): >>> i = [] >>> x, c, n, t = fsli_to_xcnt(f, s, l, i) >>> x - array([1, 2, 3, 4, 5]) + array([1., 2., 3., 4., 5.]) >>> c array([0, 1, 1, 0, 0]) >>> n @@ -1476,8 +1476,12 @@ def fs_to_xrd(f, s): >>> f = [1, 4, 5] >>> s = [2, 3] >>> x, r, d = fs_to_xrd(f, s) - >>> x, r, d - (array([1, 2, 3, 4, 5]), array([5, 4, 3, 2, 1]), array([1, 0, 0, 1, 1])) + >>> x + array([1., 2., 3., 4., 5.]) + >>> r + array([5, 4, 3, 2, 1]) + >>> d + array([1, 0, 0, 1, 1]) """ x, c, n, _ = fs_to_xcnt(f, s) return xcnt_to_xrd(x, c, n)