diff --git a/.github/CODEOWNERS b/.github/CODEOWNERS new file mode 100644 index 0000000..4e0fc6b --- /dev/null +++ b/.github/CODEOWNERS @@ -0,0 +1,2 @@ +# Every file is owned by the maintainer; pull requests need their review. +* @genedan diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index 7186afb..e6b84fe 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -10,6 +10,34 @@ permissions: contents: read jobs: + lint: + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v4 + - uses: actions/setup-python@v5 + with: + python-version: "3.12" + - name: Ruff + run: | + pip install ruff + ruff check . + ruff format --check . + + typecheck: + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v4 + - uses: actions/setup-python@v5 + with: + python-version: "3.12" + cache: pip + - name: Install + run: pip install ".[dev]" + - name: Pyright + run: pyright + - name: Type completeness + run: pyright --verifytypes fcaz --ignoreexternal + test: runs-on: ubuntu-latest strategy: diff --git a/.gitignore b/.gitignore index 70a11a4..22e853a 100644 --- a/.gitignore +++ b/.gitignore @@ -1,9 +1,10 @@ .idea/ -.venv/docs/_build/ -docs/api/generated/ +.venv/ +build/ +dist/ docs/_build/ +docs/api/generated/ *.egg-info/ __pycache__/ coverage.xml .coverage - diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml new file mode 100644 index 0000000..1fd059f --- /dev/null +++ b/.pre-commit-config.yaml @@ -0,0 +1,27 @@ +default_install_hook_types: + - pre-commit + - pre-push + +repos: + - repo: https://github.com/astral-sh/ruff-pre-commit + rev: v0.17.1 + hooks: + - id: ruff-check + - id: ruff-format + + - repo: local + hooks: + - id: pyright + name: pyright + entry: pyright + language: system + types: [python] + pass_filenames: false + stages: [pre-push] + - id: pytest + name: pytest + entry: python -m pytest -q --doctest-modules fcaz tests + language: system + types: [python] + pass_filenames: false + stages: [pre-push] diff --git a/docs/_static/.gitkeep b/docs/_static/.gitkeep new file mode 100644 index 0000000..e69de29 diff --git a/docs/api/index.rst b/docs/api/index.rst index bd22a57..91f7cf0 100644 --- a/docs/api/index.rst +++ b/docs/api/index.rst @@ -5,3 +5,4 @@ API reference :maxdepth: 1 stats + typing diff --git a/docs/api/typing.rst b/docs/api/typing.rst new file mode 100644 index 0000000..d3ad5e8 --- /dev/null +++ b/docs/api/typing.rst @@ -0,0 +1,13 @@ +fcaz.typing +=========== + +.. automodule:: fcaz.typing + +.. autosummary:: + :toctree: generated + + ArrayLikeFloat + FloatOrArray + SupportsPdfSf + SupportsLogSf + FrozenDistribution diff --git a/docs/conf.py b/docs/conf.py index f46dca2..300fa3a 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -1,7 +1,10 @@ import os import sys -sys.path.insert(0, os.path.abspath('..')) +sys.path.insert( + 0, + os.path.abspath('..'), +) from fcaz.constants import BUILD_VERSION diff --git a/fcaz/__init__.py b/fcaz/__init__.py index e69de29..e3d0d80 100644 --- a/fcaz/__init__.py +++ b/fcaz/__init__.py @@ -0,0 +1,7 @@ +"""Property and casualty actuarial modeling.""" + +from fcaz import stats + +__all__ = [ + 'stats', +] diff --git a/fcaz/constants.py b/fcaz/constants.py index aae1841..ac1a688 100644 --- a/fcaz/constants.py +++ b/fcaz/constants.py @@ -1,3 +1,5 @@ -# The build version is referenced by the setup file and documentation. +"""Package-wide constants.""" + +# The build version is referenced by pyproject.toml and the documentation. BUILD_VERSION = '0.0.0' diff --git a/fcaz/py.typed b/fcaz/py.typed new file mode 100644 index 0000000..e69de29 diff --git a/fcaz/stats.py b/fcaz/stats.py new file mode 100644 index 0000000..3ec7320 --- /dev/null +++ b/fcaz/stats.py @@ -0,0 +1,365 @@ +""" +Distribution tools that extend :mod:`scipy.stats` with actuarial transformations. + +Everything here follows the scipy conventions: ``invexpon(scale=theta)`` is used like +``scipy.stats.expon(scale=theta)``, and :class:`Mixture` returns an object with the +frozen-distribution methods (``pdf``, ``cdf``, ``sf``, ``ppf``, ``rvs``, ``mean``) plus +``hf`` and ``chf``. +""" + +import math + +from collections.abc import Sequence +from typing import Any + +import numpy as np + +from numpy.typing import NDArray +from scipy import optimize +from scipy.stats import rv_continuous + +from fcaz.typing import ( + ArrayLikeFloat, + FloatOrArray, + FrozenDistribution, + SupportsLogSf, + SupportsPdfSf, +) + +__all__ = [ + 'Mixture', + 'chf', + 'hf', + 'invexpon', +] + + +def hf( + dist: SupportsPdfSf, + x: ArrayLikeFloat, +) -> FloatOrArray: + """ + Hazard function of a frozen distribution (i.e., one with parameters already fixed). + + Parameters + ---------- + dist : frozen distribution + Any object with ``pdf`` and ``sf`` methods, e.g. ``scipy.stats.expon(scale=8)``. + x : array_like + Points at which to evaluate the hazard. + + Returns + ------- + ndarray + :math:`h(x) = f(x) / S(x)`. + + Examples + -------- + >>> from scipy import stats + >>> from fcaz.stats import hf + >>> round(float(hf(stats.expon(scale=8), 2.5)), 4) + 0.125 + """ + return np.asarray(dist.pdf(x) / dist.sf(x))[()] + + +def chf( + dist: SupportsLogSf, + x: ArrayLikeFloat, +) -> FloatOrArray: + r""" + Cumulative hazard function of a frozen distribution (i.e., one with parameters already fixed). + + Parameters + ---------- + dist : frozen distribution + Any object with a ``logsf`` method, e.g. ``scipy.stats.expon(scale=8)``. + x : array_like + Points at which to evaluate the cumulative hazard. + + Returns + ------- + ndarray + :math:`H(x) = -\\ln S(x)`. + + Examples + -------- + >>> from scipy import stats + >>> from fcaz.stats import chf + >>> round(float(chf(stats.expon(scale=8), 2.5)), 4) + 0.3125 + """ + return np.asarray(-dist.logsf(x))[()] + + +# scipy convention +# noinspection PyPep8Naming +class _invexpon_gen(rv_continuous): + r""" + Inverse exponential continuous random variable. + + %(before_notes)s + + Notes + ----- + The cumulative distribution function is + + .. math:: + + F(x) = e^{-\theta / x}, \quad x > 0, + + where :math:`\theta` is the ``scale`` parameter, matching the *Loss Models* + parameterization. This is :obj:`scipy.stats.invweibull` with ``c=1``. + + %(after_notes)s + + Examples + -------- + >>> from fcaz.stats import invexpon + >>> X = invexpon(scale=4) + >>> round(float(X.cdf(6)), 4) + 0.5134 + >>> round(float(X.ppf(0.3)), 4) + 3.3223 + """ + + def _pdf( + self, + x: Any, + *args: Any, + ) -> Any: + """Probability density function of the standard form, ``exp(-1/x) / x**2``.""" + return np.exp(-1.0 / x) / x**2 + + def _cdf( + self, + x: Any, + *args: Any, + ) -> Any: + """Cumulative distribution function of the standard form, ``exp(-1/x)``.""" + return np.exp(-1.0 / x) + + def _sf( + self, + x: Any, + *args: Any, + ) -> Any: + """Survival function of the standard form, ``1 - exp(-1/x)``.""" + return -np.expm1(-1.0 / x) + + def _ppf( + self, + q: Any, + *args: Any, + ) -> Any: + """Percent point function of the standard form, ``-1 / log(q)``.""" + return -1.0 / np.log(q) + + +# Inverse exponential distribution from Loss Models. +invexpon: rv_continuous = _invexpon_gen( + a=0.0, + name='invexpon', +) + + +class Mixture: + r""" + Finite mixture of frozen distributions (i.e., those with parameters already fixed). + + Parameters + ---------- + dists : sequence of frozen distributions + Component distributions, e.g. ``scipy.stats.expon(scale=8)``. + weights : array_like + Non-negative mixing weights, one per component, summing to 1. + + Raises + ------ + ValueError + If ``dists`` and ``weights`` differ in length, or the weights are negative or do not + sum to 1. + + Notes + ----- + The survival function is :math:`S(x) = \sum_i w_i S_i(x)`; ``pdf``, ``cdf`` and + ``mean`` mix in the same way. ``hf`` and ``chf`` are derived from the mixed ``pdf`` and + ``sf``. + + Examples + -------- + A 40% exponential, 60% inverse exponential mixture: + + >>> from scipy import stats + >>> from fcaz.stats import invexpon, Mixture + >>> X = Mixture([stats.expon(scale=8), invexpon(scale=3)], [0.4, 0.6]) + >>> round(float(X.chf(2.5)), 4) + 0.3398 + """ + + def __init__( + self, + dists: Sequence[FrozenDistribution], + weights: ArrayLikeFloat, + ) -> None: + w: NDArray[np.float64] = np.asarray( + weights, + dtype=np.float64, + ) + if len(dists) != len(w): + raise ValueError('dists and weights must have the same length') + if bool(w.min() < 0) or not math.isclose( + float(w.sum()), + 1.0, + rel_tol=1e-5, + abs_tol=1e-8, + ): + raise ValueError('weights must be non-negative and sum to 1') + self.dists: list[FrozenDistribution] = list(dists) + self.weights: NDArray[np.float64] = w + + def _mix( + self, + method: str, + x: ArrayLikeFloat, + ) -> FloatOrArray: + r""" + Weighted sum of a method evaluated on every component. + + Parameters + ---------- + method : str + Name of the component method to call, e.g. ``'pdf'``, ``'cdf'`` or ``'sf'``. + x : array_like + Points at which to evaluate the method. + + Returns + ------- + float or ndarray + :math:`\sum_i w_i \, g_i(x)`, where :math:`g_i` is ``method`` of component + :math:`i`. Valid for any quantity that mixes linearly, which excludes ``ppf`` + and ``logsf``. + """ + xs: NDArray[np.float64] = np.asarray( + x, + dtype=np.float64, + ) + total: NDArray[np.float64] = np.zeros_like(xs) + for w, d in zip( + self.weights, + self.dists, + strict=True, + ): + total = total + w * np.asarray( + getattr( + d, + method, + )(xs), + dtype=np.float64, + ) + return total[()] + + def pdf(self, x: ArrayLikeFloat) -> FloatOrArray: + """Probability density function.""" + return self._mix( + 'pdf', + x, + ) + + def cdf(self, x: ArrayLikeFloat) -> FloatOrArray: + """Cumulative distribution function.""" + return self._mix( + 'cdf', + x, + ) + + def sf(self, x: ArrayLikeFloat) -> FloatOrArray: + """Survival function, 1 - cdf.""" + return self._mix( + 'sf', + x, + ) + + def logsf(self, x: ArrayLikeFloat) -> FloatOrArray: + """Log of the survival function.""" + return np.log(self.sf(x)) + + def hf(self, x: ArrayLikeFloat) -> FloatOrArray: + """Hazard function, pdf / sf.""" + return hf( + self, + x, + ) + + def chf(self, x: ArrayLikeFloat) -> FloatOrArray: + """Cumulative hazard function, -log sf.""" + return chf( + self, + x, + ) + + def mean(self) -> float: + """Mean of the mixture, the weighted mean of the component means.""" + return float( + sum( + w * d.mean() + for w, d in zip( + self.weights, + self.dists, + strict=True, + ) + ) + ) + + def ppf(self, q: ArrayLikeFloat) -> FloatOrArray: + """Percent point function (inverse cdf), solved numerically.""" + + def solve(p: float) -> float: + lo = min(d.ppf(p) for d in self.dists) + hi = max(d.ppf(p) for d in self.dists) + return ( + lo + if lo == hi + else optimize.brentq( + lambda t: float(self.cdf(t)) - p, + lo, + hi, + ) + ) + + return np.vectorize( + solve, + otypes=[float], + )(q) + + def rvs( + self, + size: int | tuple[int, ...] | None = None, + random_state: Any = None, + ) -> FloatOrArray: + """Random variates; ``size`` and ``random_state`` as in scipy.""" + rng = np.random.default_rng(random_state) + n = 1 if size is None else int(np.prod(size)) + edges = self.weights.cumsum() + which = np.minimum( + (rng.random(n)[:, None] >= edges).sum(axis=1), + len(self.dists) - 1, + ) + out: NDArray[np.float64] = np.empty(n) + for i, d in enumerate(self.dists): + m = which == i + out[m] = d.rvs( + size=int(m.sum()), + random_state=rng, + ) + if size is None: + return out[0] + shape: tuple[int, ...] = ( + (size,) + if isinstance( + size, + int, + ) + else size + ) + return out.reshape(shape) # pyright: ignore[reportUnknownMemberType] diff --git a/fcaz/typing.py b/fcaz/typing.py new file mode 100644 index 0000000..ce200b3 --- /dev/null +++ b/fcaz/typing.py @@ -0,0 +1,78 @@ +""" +Type aliases and protocols shared across fcaz. + +``FrozenDistribution`` and the smaller ``Supports*`` protocols describe the methods fcaz needs +from a frozen distribution, so any frozen ``scipy.stats`` distribution or fcaz object can be +passed where one is expected. +""" + +from collections.abc import Sequence +from typing import ( + Any, + Protocol, + TypeAlias, +) + +import numpy as np + +from numpy.typing import NDArray + +__all__ = [ + 'ArrayLikeFloat', + 'FloatOrArray', + 'FrozenDistribution', + 'SupportsLogSf', + 'SupportsPdfSf', +] + +ArrayLikeFloat: TypeAlias = ( + float + | np.floating[Any] + | NDArray[np.number[Any]] + | Sequence[float | np.floating[Any] | NDArray[np.number[Any]]] +) +FloatOrArray: TypeAlias = np.float64 | NDArray[np.float64] + + +class SupportsPdfSf(Protocol): + """Any object with ``pdf`` and ``sf`` methods, such as a frozen scipy distribution.""" + + def pdf(self, x: Any) -> Any: + """Probability density function.""" + ... + + def sf(self, x: Any) -> Any: + """Survival function.""" + ... + + +class SupportsLogSf(Protocol): + """Any object with a ``logsf`` method, such as a frozen scipy distribution.""" + + def logsf(self, x: Any) -> Any: + """Log of the survival function.""" + ... + + +class FrozenDistribution(SupportsPdfSf, SupportsLogSf, Protocol): + """The frozen-distribution interface that Mixture needs from its components.""" + + def cdf(self, x: Any) -> Any: + """Cumulative distribution function.""" + ... + + def ppf(self, q: Any) -> Any: + """Percent point function.""" + ... + + def mean(self) -> Any: + """Mean.""" + ... + + def rvs( + self, + size: Any = None, + random_state: Any = None, + ) -> Any: + """Random variates.""" + ... diff --git a/pyproject.toml b/pyproject.toml index f0e43c4..658a856 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -19,13 +19,20 @@ classifiers = [ [project.optional-dependencies] dev = [ "build", + "pre-commit", + "pyright", "pytest", + "ruff", + "scipy-stubs", "twine" ] docs = [ "sphinx", "pydata-sphinx-theme" ] +all = [ + "fcaz[dev, docs]" +] [project.urls] Homepage = "https://github.com/casact/fcaz" @@ -43,3 +50,59 @@ branch = true [tool.coverage.report] show_missing = true + +[tool.ruff] +line-length = 100 +target-version = "py310" +extend-exclude = ["docs/_build"] + +[tool.ruff.lint] +select = [ + "E", + "F", + "W", + "I", + "B", + "UP", + "NPY", + "SIM", + "RUF", + "D", + "D213", +] +ignore = ["D212"] + +[tool.ruff.lint.pydocstyle] +convention = "numpy" + +[tool.ruff.lint.isort] +known-first-party = ["fcaz"] +force-sort-within-sections = false +order-by-type = true +combine-as-imports = false +lines-after-imports = -1 +lines-between-types = 1 +split-on-trailing-comma = true + +[tool.ruff.lint.per-file-ignores] +"tests/*" = ["D"] +"docs/conf.py" = ["D"] + +[tool.ruff.format] +quote-style = "single" +docstring-code-format = true + +[tool.setuptools.package-data] +fcaz = ["py.typed"] + +[tool.pyright] +include = ["fcaz", "tests"] +exclude = ["docs/_build"] +pythonVersion = "3.10" +strict = ["fcaz"] + +[[tool.pyright.executionEnvironments]] +root = "tests" +extraPaths = ["."] +reportMissingParameterType = "error" +reportMissingTypeArgument = "error" diff --git a/tests/test_stats.py b/tests/test_stats.py new file mode 100644 index 0000000..e26820c --- /dev/null +++ b/tests/test_stats.py @@ -0,0 +1,205 @@ +from typing import TypedDict + +import numpy as np +import pytest + +from numpy.typing import ( + ArrayLike, + NDArray, +) +from scipy import stats + +from fcaz.stats import ( + Mixture, + chf, + hf, + invexpon, +) +from fcaz.typing import FrozenDistribution + +SEEDS = range(20) + + +class Draw(TypedDict): + """Random inputs shared by the tests.""" + + theta1: float | int + theta2: float | int + w: float | int + x: NDArray[np.float64] + q: NDArray[np.float64] + + +def close( + a: ArrayLike, + b: ArrayLike, +) -> bool: + """Whether every element of a is numerically close to b.""" + return bool(np.allclose(a, b)) + + +def draw(seed: int) -> Draw: + """Random scales, mixing weight, evaluation points and probabilities for a seed.""" + rng = np.random.default_rng(seed) + return { + 'theta1': rng.uniform( + 0.5, + 20, + ), + 'theta2': rng.uniform( + 0.5, + 20, + ), + 'w': rng.uniform( + 0.05, + 0.95, + ), + 'x': rng.uniform( + 0.1, + 30, + size=5, + ), + 'q': rng.uniform( + 0.01, + 0.99, + size=5, + ), + } + + +@pytest.mark.parametrize( + 'seed', + SEEDS, +) +def test_invexpon_matches_invweibull_and_closed_form(seed: int) -> None: + """invexpon agrees with its closed-form cdf and with scipy's invweibull at c=1.""" + p = draw(seed) + th, x, q = p['theta1'], p['x'], p['q'] + d: FrozenDistribution = invexpon(scale=th) + ref: FrozenDistribution = stats.invweibull( + c=1, + scale=th, + ) + assert close( + d.cdf(x), + np.exp(-th / x), + ) + assert close( + d.cdf(x), + ref.cdf(x), + ) + assert close( + d.pdf(x), + ref.pdf(x), + ) + assert close( + d.sf(x), + ref.sf(x), + ) + assert close( + d.ppf(q), + ref.ppf(q), + ) + + +@pytest.mark.parametrize( + 'seed', + SEEDS, +) +def test_mixture_sf_and_chf_closed_form(seed: int) -> None: + """Mixture sf, cdf, chf and hf match closed forms for an exponential and invexpon mix.""" + p = draw(seed) + t1, t2, w, x = p['theta1'], p['theta2'], p['w'], p['x'] + m = Mixture( + [stats.expon(scale=t1), invexpon(scale=t2)], + [w, 1 - w], + ) + s = w * np.exp(-x / t1) + (1 - w) * (1 - np.exp(-t2 / x)) + f = w * np.exp(-x / t1) / t1 + (1 - w) * np.exp(-t2 / x) * t2 / x**2 + assert close( + m.sf(x), + s, + ) + assert close( + m.cdf(x), + 1 - s, + ) + assert close( + m.chf(x), + -np.log(s), + ) + assert close( + m.hf(x), + f / s, + ) + + +@pytest.mark.parametrize( + 'seed', + SEEDS, +) +def test_hf_chf_helpers_on_exponential(seed: int) -> None: + """hf and chf of an exponential give the hazard 1/theta and cumulative hazard x/theta.""" + p = draw(seed) + th, x = p['theta1'], p['x'] + e: FrozenDistribution = stats.expon(scale=th) + assert close( + chf( + e, + x, + ), + x / th, + ) + assert close( + hf( + e, + x, + ), + 1 / th, + ) + + +@pytest.mark.parametrize( + 'seed', + SEEDS, +) +def test_mixture_ppf_roundtrip_and_mean(seed: int) -> None: + """Mixture ppf inverts cdf, and the mean and sample mean match the weighted component means.""" + p = draw(seed) + t1, t2, w, q = p['theta1'], p['theta2'], p['w'], p['q'] + m = Mixture( + [stats.expon(scale=t1), invexpon(scale=t2)], + [w, 1 - w], + ) + assert close( + m.cdf(m.ppf(q)), + q, + ) + m2 = Mixture( + [stats.expon(scale=t1), stats.expon(scale=t2)], + [w, 1 - w], + ) + assert close( + m2.mean(), + w * t1 + (1 - w) * t2, + ) + r = m2.rvs( + size=100_000, + random_state=seed, + ) + assert abs(r.mean() - m2.mean()) < 0.05 * m2.mean() + + +def test_mixture_rejects_bad_weights() -> None: + """Mixture raises ValueError for weights that do not sum to 1 or have the wrong length.""" + e: FrozenDistribution = stats.expon() + with pytest.raises(ValueError): + Mixture( + [e, e], + [0.5, 0.6], + ) + with pytest.raises(ValueError): + Mixture( + [e, e], + [1.0], + )