Skip to content

lossmodels

Severity and frequency distributions, maximum-likelihood fitting, goodness-of-fit, and collective-risk (compound) models — the loss-distribution toolkit familiar from the actuarial FAM/STAM syllabus, as composable objects.

Quickstart

import lossmodels as lm

# a severity distribution
sev = lm.Lognormal(mu=8.0, sigma=1.5)
sev.quantile(0.99)                     # 99th-percentile claim
sev.limited_expected_value(250_000)    # E[min(X, 250k)] -- a capped/limited loss

# a collective-risk (compound) model: random number of claims, random sizes
model = lm.CollectiveRiskModel(
    frequency=lm.Poisson(lam=3.0),
    severity=lm.Lognormal(mu=8.0, sigma=1.5),
)
model.mean()                           # E[aggregate]
model.var(0.99)                        # 99% Value-at-Risk of the aggregate
model.tvar(0.99)                       # 99% Tail VaR

Fitting

import numpy as np
import lossmodels as lm

data = lm.Lognormal(mu=8.0, sigma=1.2).sample(5_000, seed=0)
fit = lm.fit_lognormal(data)           # MLE
lm.fit_best_severity(data)             # compare candidates by AIC/BIC

API reference

lossmodels

lossmodels: actuarial loss distributions, aggregate modeling, and fit diagnostics.

The full public API is surfaced here for convenience, so common entry points are available directly as lossmodels.Lognormal, lossmodels.fit_best_severity, lossmodels.goodness_of_fit, etc. Submodule imports (from lossmodels.severity import Lognormal) continue to work unchanged.

SeverityModel

Bases: ABC

Base class for severity (loss size) distributions.

All severity models must implement:

  • sample
  • mean
  • variance
Source code in lossmodels/severity/base.py 
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
class SeverityModel(ABC):
    """
    Base class for severity (loss size) distributions.

    All severity models must implement:

    - sample
    - mean
    - variance
    """

    @abstractmethod
    def sample(self, size: int = 1) -> np.ndarray:
        """Generate random loss samples"""
        pass

    @abstractmethod
    def mean(self) -> float:
        """Expected loss"""
        pass

    @abstractmethod
    def variance(self) -> float:
        """Variance of loss"""
        pass

    def std(self) -> float:
        return np.sqrt(self.variance())

    # --- Quantile / inverse CDF ---
    def quantile(self, p):
        """Inverse CDF (quantile / Value-at-Risk).

        Numerically inverts ``cdf``. Subclasses with a closed-form inverse
        override this for speed and accuracy. Returns a Python ``float`` for a
        scalar ``p`` and a ``numpy.ndarray`` for array-like ``p``.
        """
        arr = np.asarray(p, dtype=float)
        out = np.array([self._invert_cdf(float(pi)) for pi in np.atleast_1d(arr)], dtype=float)
        return float(out[0]) if arr.ndim == 0 else out.reshape(arr.shape)

    def ppf(self, p):
        """Alias for :meth:`quantile` (inverse CDF)."""
        return self.quantile(p)

    def _invert_cdf(self, p: float) -> float:
        if not 0.0 <= p <= 1.0:
            raise ValueError("p must be in [0, 1].")
        if p == 0.0:
            return 0.0
        if p == 1.0:
            return float("inf")
        if not hasattr(self, "cdf"):
            raise TypeError("Severity model must implement cdf(x) to invert.")
        hi = 1.0
        while self.cdf(hi) < p:
            hi *= 2.0
            if hi > 1e15:
                break
        return float(brentq(lambda x: self.cdf(x) - p, 0.0, hi, maxiter=200))

    # --- Actuarial-specific methods ---
    def limited_expected_value(self, d: float, n_sim: int = 100_000) -> float:
        """
        E[min(X, d)] computed deterministically from the survival function.

        For a nonnegative loss random variable X,

            E[min(X, d)] = integral_0^d S_X(x) dx

        where S_X(x) = 1 - F_X(x).

        The ``n_sim`` argument is retained for backward compatibility but is no
        longer used.
        """
        del n_sim  # backward-compatibility placeholder

        if d < 0:
            raise ValueError("d must be nonnegative.")
        if d == 0:
            return 0.0
        if not hasattr(self, "cdf"):
            raise TypeError("Severity model must implement cdf(x).")

        value, _ = quad(lambda x: 1.0 - self.cdf(x), 0.0, d, limit=200)
        return float(value)

    def excess_loss(self, d: float, n_sim: int = 100_000) -> float:
        """
        E[(X - d)+] computed deterministically from the survival function.

        For a nonnegative loss random variable X,

            E[(X - d)+] = integral_d^infinity S_X(x) dx

        where S_X(x) = 1 - F_X(x).

        This remains well defined even when E[X] does not exist, provided the
        tail integral from d to infinity converges. The ``n_sim`` argument is
        retained for backward compatibility but is no longer used.
        """
        del n_sim  # backward-compatibility placeholder

        if d < 0:
            raise ValueError("d must be nonnegative.")
        if not hasattr(self, "cdf"):
            raise TypeError("Severity model must implement cdf(x).")

        value, _ = quad(lambda x: 1.0 - self.cdf(x), d, np.inf, limit=200)
        return float(value)

    def __repr__(self):
        return f"{self.__class__.__name__}()"
sample abstractmethod 
sample(size: int = 1) -> np.ndarray

Generate random loss samples

Source code in lossmodels/severity/base.py 
19
20
21
22
@abstractmethod
def sample(self, size: int = 1) -> np.ndarray:
    """Generate random loss samples"""
    pass
mean abstractmethod 
mean() -> float

Expected loss

Source code in lossmodels/severity/base.py 
24
25
26
27
@abstractmethod
def mean(self) -> float:
    """Expected loss"""
    pass
variance abstractmethod 
variance() -> float

Variance of loss

Source code in lossmodels/severity/base.py 
29
30
31
32
@abstractmethod
def variance(self) -> float:
    """Variance of loss"""
    pass
quantile 
quantile(p)

Inverse CDF (quantile / Value-at-Risk).

Numerically inverts cdf. Subclasses with a closed-form inverse override this for speed and accuracy. Returns a Python float for a scalar p and a numpy.ndarray for array-like p.

Source code in lossmodels/severity/base.py 
38
39
40
41
42
43
44
45
46
47
def quantile(self, p):
    """Inverse CDF (quantile / Value-at-Risk).

    Numerically inverts ``cdf``. Subclasses with a closed-form inverse
    override this for speed and accuracy. Returns a Python ``float`` for a
    scalar ``p`` and a ``numpy.ndarray`` for array-like ``p``.
    """
    arr = np.asarray(p, dtype=float)
    out = np.array([self._invert_cdf(float(pi)) for pi in np.atleast_1d(arr)], dtype=float)
    return float(out[0]) if arr.ndim == 0 else out.reshape(arr.shape)
ppf 
ppf(p)

Alias for :meth:quantile (inverse CDF).

Source code in lossmodels/severity/base.py 
49
50
51
def ppf(self, p):
    """Alias for :meth:`quantile` (inverse CDF)."""
    return self.quantile(p)
limited_expected_value 
limited_expected_value(
    d: float, n_sim: int = 100000
) -> float

E[min(X, d)] computed deterministically from the survival function.

For a nonnegative loss random variable X,

E[min(X, d)] = integral_0^d S_X(x) dx

where S_X(x) = 1 - F_X(x).

The n_sim argument is retained for backward compatibility but is no longer used.

Source code in lossmodels/severity/base.py 
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
def limited_expected_value(self, d: float, n_sim: int = 100_000) -> float:
    """
    E[min(X, d)] computed deterministically from the survival function.

    For a nonnegative loss random variable X,

        E[min(X, d)] = integral_0^d S_X(x) dx

    where S_X(x) = 1 - F_X(x).

    The ``n_sim`` argument is retained for backward compatibility but is no
    longer used.
    """
    del n_sim  # backward-compatibility placeholder

    if d < 0:
        raise ValueError("d must be nonnegative.")
    if d == 0:
        return 0.0
    if not hasattr(self, "cdf"):
        raise TypeError("Severity model must implement cdf(x).")

    value, _ = quad(lambda x: 1.0 - self.cdf(x), 0.0, d, limit=200)
    return float(value)
excess_loss 
excess_loss(d: float, n_sim: int = 100000) -> float

E[(X - d)+] computed deterministically from the survival function.

For a nonnegative loss random variable X,

E[(X - d)+] = integral_d^infinity S_X(x) dx

where S_X(x) = 1 - F_X(x).

This remains well defined even when E[X] does not exist, provided the tail integral from d to infinity converges. The n_sim argument is retained for backward compatibility but is no longer used.

Source code in lossmodels/severity/base.py 
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
def excess_loss(self, d: float, n_sim: int = 100_000) -> float:
    """
    E[(X - d)+] computed deterministically from the survival function.

    For a nonnegative loss random variable X,

        E[(X - d)+] = integral_d^infinity S_X(x) dx

    where S_X(x) = 1 - F_X(x).

    This remains well defined even when E[X] does not exist, provided the
    tail integral from d to infinity converges. The ``n_sim`` argument is
    retained for backward compatibility but is no longer used.
    """
    del n_sim  # backward-compatibility placeholder

    if d < 0:
        raise ValueError("d must be nonnegative.")
    if not hasattr(self, "cdf"):
        raise TypeError("Severity model must implement cdf(x).")

    value, _ = quad(lambda x: 1.0 - self.cdf(x), d, np.inf, limit=200)
    return float(value)

Exponential

Bases: SeverityModel

Exponential severity model.

Parameterization

X ~ Exponential(rate)

Support: x >= 0

Mean = 1 / rate Variance = 1 / rate^2

Parameters

rate : float Rate parameter (lambda), with rate > 0.

Source code in lossmodels/severity/exponential.py 
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
class Exponential(SeverityModel):
    """
    Exponential severity model.

    Parameterization
    ----------------
    X ~ Exponential(rate)

    Support: x >= 0

    Mean = 1 / rate
    Variance = 1 / rate^2

    Parameters
    ----------
    rate : float
        Rate parameter (lambda), with rate > 0.
    """

    def __init__(self, rate: float):
        if rate <= 0:
            raise ValueError("rate must be positive.")

        self.rate = rate
        self.scale = 1.0 / rate  # SciPy uses scale = 1 / rate

    def sample(self, size: int = 1) -> np.ndarray:
        """
        Generate random samples.
        """
        if size <= 0:
            raise ValueError("size must be positive.")

        return np.random.exponential(scale=self.scale, size=size)

    def mean(self) -> float:
        return 1.0 / self.rate

    def variance(self) -> float:
        return 1.0 / (self.rate ** 2)

    def pdf(self, x):
        return eval_dist(lambda v: expon.pdf(v, scale=self.scale), x)

    def cdf(self, x):
        return eval_dist(lambda v: expon.cdf(v, scale=self.scale), x)

    def quantile(self, p):
        return eval_dist(lambda v: expon.ppf(v, scale=self.scale), p)

    def excess_loss(self, d: float) -> float:
        """
        E[(X - d)+] = exp(-rate * d) / rate
        """
        if d < 0:
            raise ValueError("d must be nonnegative.")

        return float(np.exp(-self.rate * d) / self.rate)

    def limited_expected_value(self, d: float) -> float:
        """
        E[min(X, d)] = (1 - exp(-rate * d)) / rate
        """
        if d < 0:
            raise ValueError("d must be nonnegative.")

        return float((1.0 - np.exp(-self.rate * d)) / self.rate)

    def __repr__(self) -> str:
        return f"Exponential(rate={self.rate})"
sample 
sample(size: int = 1) -> np.ndarray

Generate random samples.

Source code in lossmodels/severity/exponential.py 
34
35
36
37
38
39
40
41
def sample(self, size: int = 1) -> np.ndarray:
    """
    Generate random samples.
    """
    if size <= 0:
        raise ValueError("size must be positive.")

    return np.random.exponential(scale=self.scale, size=size)
excess_loss 
excess_loss(d: float) -> float

E[(X - d)+] = exp(-rate * d) / rate

Source code in lossmodels/severity/exponential.py 
58
59
60
61
62
63
64
65
def excess_loss(self, d: float) -> float:
    """
    E[(X - d)+] = exp(-rate * d) / rate
    """
    if d < 0:
        raise ValueError("d must be nonnegative.")

    return float(np.exp(-self.rate * d) / self.rate)
limited_expected_value 
limited_expected_value(d: float) -> float

E[min(X, d)] = (1 - exp(-rate * d)) / rate

Source code in lossmodels/severity/exponential.py 
67
68
69
70
71
72
73
74
def limited_expected_value(self, d: float) -> float:
    """
    E[min(X, d)] = (1 - exp(-rate * d)) / rate
    """
    if d < 0:
        raise ValueError("d must be nonnegative.")

    return float((1.0 - np.exp(-self.rate * d)) / self.rate)

Gamma

Bases: SeverityModel

Gamma severity model.

Parameterization

X ~ Gamma(shape=alpha, scale=theta) Support: x > 0

Parameters

alpha : float Shape parameter, with alpha > 0. theta : float Scale parameter, with theta > 0.

Source code in lossmodels/severity/gamma.py 
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
class Gamma(SeverityModel):
    """
    Gamma severity model.

    Parameterization
    ----------------
    X ~ Gamma(shape=alpha, scale=theta)
    Support: x > 0

    Parameters
    ----------
    alpha : float
        Shape parameter, with alpha > 0.
    theta : float
        Scale parameter, with theta > 0.
    """

    def __init__(self, alpha: float, theta: float):
        if alpha <= 0:
            raise ValueError("alpha must be positive.")
        if theta <= 0:
            raise ValueError("theta must be positive.")

        self.alpha = alpha
        self.theta = theta

    def sample(self, size: int = 1) -> np.ndarray:
        if size <= 0:
            raise ValueError("size must be positive.")

        return np.random.gamma(shape=self.alpha, scale=self.theta, size=size)

    def mean(self) -> float:
        return self.alpha * self.theta

    def variance(self) -> float:
        return self.alpha * (self.theta ** 2)

    def pdf(self, x):
        return eval_dist(lambda v: gamma.pdf(v, a=self.alpha, scale=self.theta), x)

    def cdf(self, x):
        return eval_dist(lambda v: gamma.cdf(v, a=self.alpha, scale=self.theta), x)

    def quantile(self, p):
        return eval_dist(lambda v: gamma.ppf(v, a=self.alpha, scale=self.theta), p)

    def __repr__(self) -> str:
        return f"Gamma(alpha={self.alpha}, theta={self.theta})"

Lognormal

Bases: SeverityModel

Lognormal severity model.

Parameterization

If Y = log(X) ~ Normal(mu, sigma^2), then X is Lognormal(mu, sigma). Support: x > 0

Parameters

mu : float Mean of log(X). sigma : float Standard deviation of log(X), with sigma > 0.

Source code in lossmodels/severity/lognormal.py 
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
class Lognormal(SeverityModel):
    """
    Lognormal severity model.

    Parameterization
    ----------------
    If Y = log(X) ~ Normal(mu, sigma^2), then X is Lognormal(mu, sigma).
    Support: x > 0

    Parameters
    ----------
    mu : float
        Mean of log(X).
    sigma : float
        Standard deviation of log(X), with sigma > 0.
    """

    def __init__(self, mu: float, sigma: float):
        if sigma <= 0:
            raise ValueError("sigma must be positive.")

        self.mu = mu
        self.sigma = sigma

    def sample(self, size: int = 1) -> np.ndarray:
        if size <= 0:
            raise ValueError("size must be positive.")

        return np.random.lognormal(mean=self.mu, sigma=self.sigma, size=size)

    def mean(self) -> float:
        return float(np.exp(self.mu + 0.5 * self.sigma ** 2))

    def variance(self) -> float:
        sigma2 = self.sigma ** 2
        return float((np.exp(sigma2) - 1) * np.exp(2 * self.mu + sigma2))

    def pdf(self, x):
        return eval_dist(lambda v: lognorm.pdf(v, s=self.sigma, scale=np.exp(self.mu)), x)

    def cdf(self, x):
        return eval_dist(lambda v: lognorm.cdf(v, s=self.sigma, scale=np.exp(self.mu)), x)

    def quantile(self, p):
        return eval_dist(lambda v: lognorm.ppf(v, s=self.sigma, scale=np.exp(self.mu)), p)

    def __repr__(self) -> str:
        return f"Lognormal(mu={self.mu}, sigma={self.sigma})"

Pareto

Bases: SeverityModel

Pareto Type I severity model.

Parameterization

X ~ Pareto(alpha, theta) Support: x >= theta

Density

f(x) = alpha * theta^alpha / x^(alpha + 1), x >= theta

Parameters

alpha : float Shape parameter, with alpha > 0. theta : float Scale (minimum) parameter, with theta > 0.

Source code in lossmodels/severity/pareto.py 
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
class Pareto(SeverityModel):
    """
    Pareto Type I severity model.

    Parameterization
    ----------------
    X ~ Pareto(alpha, theta)
    Support: x >= theta

    Density:
        f(x) = alpha * theta^alpha / x^(alpha + 1),  x >= theta

    Parameters
    ----------
    alpha : float
        Shape parameter, with alpha > 0.
    theta : float
        Scale (minimum) parameter, with theta > 0.
    """

    def __init__(self, alpha: float, theta: float):
        if alpha <= 0:
            raise ValueError("alpha must be positive.")
        if theta <= 0:
            raise ValueError("theta must be positive.")

        self.alpha = alpha
        self.theta = theta

    def sample(self, size: int = 1) -> np.ndarray:
        if size <= 0:
            raise ValueError("size must be positive.")

        return pareto.rvs(b=self.alpha, scale=self.theta, size=size)

    def mean(self) -> float:
        if self.alpha <= 1:
            raise ValueError("Mean does not exist for alpha <= 1.")
        return self.alpha * self.theta / (self.alpha - 1)

    def variance(self) -> float:
        if self.alpha <= 2:
            raise ValueError("Variance does not exist for alpha <= 2.")
        return (self.alpha * self.theta ** 2) / ((self.alpha - 1) ** 2 * (self.alpha - 2))

    def pdf(self, x):
        return eval_dist(lambda v: pareto.pdf(v, b=self.alpha, scale=self.theta), x)

    def cdf(self, x):
        return eval_dist(lambda v: pareto.cdf(v, b=self.alpha, scale=self.theta), x)

    def quantile(self, p):
        return eval_dist(lambda v: pareto.ppf(v, b=self.alpha, scale=self.theta), p)

    def __repr__(self) -> str:
        return f"Pareto(alpha={self.alpha}, theta={self.theta})"

SplicedSeverity

Bases: SeverityModel

Two-piece spliced severity: a body below a threshold, a tail above it.

The density is

f(x) = w * f_body(x) / F_body(u),    0 < x <= u
f(x) = (1 - w) * f_tail(x),          x > u

where u is the threshold and w = P(X <= u) is the body mass. The body is renormalized onto (0, u] and the tail must be a distribution supported on [u, inf) with F_tail(u) = 0 (e.g. a GPD with location u from extremeloss, or a :class:~lossmodels.Pareto with theta = u). The CDF is

F(x) = w * F_body(x) / F_body(u),    x <= u
F(x) = w + (1 - w) * F_tail(x),      x > u

which is continuous at u (both pieces equal w there). The density may jump at u; this is the ordinary (unconstrained) spliced model, not a smooth splice.

The result is an ordinary severity model -- it exposes sample(size) and mean(), so it drops straight back into risksim and extremeloss as a tail-corrected severity.

Parameters

body : severity model Body distribution. Must implement cdf, pdf, quantile and limited_expected_value (every lossmodels severity qualifies). tail : tail distribution on [u, inf) Must implement cdf (with cdf(u) == 0), pdf and sample(size). quantile enables an exact spliced quantile; mean / variance enable the spliced moments. threshold : float Split point u > 0. weight : float Body mass w = P(X <= u) in (0, 1).

Source code in lossmodels/severity/spliced.py 
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
class SplicedSeverity(SeverityModel):
    """Two-piece spliced severity: a body below a threshold, a tail above it.

    The density is

        f(x) = w * f_body(x) / F_body(u),    0 < x <= u
        f(x) = (1 - w) * f_tail(x),          x > u

    where ``u`` is the threshold and ``w = P(X <= u)`` is the body mass. The
    body is renormalized onto ``(0, u]`` and the tail must be a distribution
    supported on ``[u, inf)`` with ``F_tail(u) = 0`` (e.g. a GPD with location
    ``u`` from ``extremeloss``, or a :class:`~lossmodels.Pareto` with
    ``theta = u``). The CDF is

        F(x) = w * F_body(x) / F_body(u),    x <= u
        F(x) = w + (1 - w) * F_tail(x),      x > u

    which is continuous at ``u`` (both pieces equal ``w`` there). The density
    may jump at ``u``; this is the ordinary (unconstrained) spliced model, not a
    smooth splice.

    The result is an ordinary severity model -- it exposes ``sample(size)`` and
    ``mean()``, so it drops straight back into ``risksim`` and ``extremeloss``
    as a tail-corrected severity.

    Parameters
    ----------
    body : severity model
        Body distribution. Must implement ``cdf``, ``pdf``, ``quantile`` and
        ``limited_expected_value`` (every ``lossmodels`` severity qualifies).
    tail : tail distribution on ``[u, inf)``
        Must implement ``cdf`` (with ``cdf(u) == 0``), ``pdf`` and
        ``sample(size)``. ``quantile`` enables an exact spliced quantile;
        ``mean`` / ``variance`` enable the spliced moments.
    threshold : float
        Split point ``u > 0``.
    weight : float
        Body mass ``w = P(X <= u)`` in ``(0, 1)``.
    """

    def __init__(self, body, tail, threshold: float, weight: float):
        if threshold <= 0:
            raise ValueError("threshold must be positive.")
        if not 0.0 < weight < 1.0:
            raise ValueError("weight must be in (0, 1).")
        for label, obj, required in (
            ("body", body, ("cdf", "pdf", "quantile")),
            ("tail", tail, ("cdf", "pdf", "sample")),
        ):
            missing = [m for m in required if not hasattr(obj, m)]
            if missing:
                raise TypeError(f"{label} is missing required method(s): {', '.join(missing)}")

        self.body = body
        self.tail = tail
        self.threshold = float(threshold)
        self.weight = float(weight)

        self._Fb_u = float(body.cdf(self.threshold))
        if self._Fb_u <= 0.0:
            raise ValueError("body must place positive probability below the threshold.")
        tail_at_u = float(tail.cdf(self.threshold))
        if abs(tail_at_u) > 1e-6:
            raise ValueError(
                "tail must be supported on [threshold, inf) with cdf(threshold)=0; "
                f"got cdf(threshold)={tail_at_u:.3g}."
            )

    # --- distribution functions ---
    def cdf(self, x):
        arr = np.asarray(x, dtype=float)
        xx = np.atleast_1d(arr)
        below = xx <= self.threshold
        fb = _to_array(self.body.cdf, xx)
        ft = _to_array(self.tail.cdf, xx)
        out = np.where(
            below,
            self.weight * fb / self._Fb_u,
            self.weight + (1.0 - self.weight) * ft,
        )
        out = np.clip(out, 0.0, 1.0)
        return float(out[0]) if arr.ndim == 0 else out.reshape(arr.shape)

    def pdf(self, x):
        arr = np.asarray(x, dtype=float)
        xx = np.atleast_1d(arr)
        below = xx <= self.threshold
        fb = _to_array(self.body.pdf, xx)
        ft = _to_array(self.tail.pdf, xx)
        out = np.where(
            below,
            self.weight * fb / self._Fb_u,
            (1.0 - self.weight) * ft,
        )
        return float(out[0]) if arr.ndim == 0 else out.reshape(arr.shape)

    def quantile(self, p):
        arr = np.asarray(p, dtype=float)
        ps = np.atleast_1d(arr)
        out = np.empty_like(ps)
        lo = ps <= self.weight
        if lo.any():
            qb = self.body.quantile(ps[lo] * self._Fb_u / self.weight)
            out[lo] = np.atleast_1d(np.asarray(qb, dtype=float))
        hi = ~lo
        if hi.any():
            pt = (ps[hi] - self.weight) / (1.0 - self.weight)
            if hasattr(self.tail, "quantile"):
                out[hi] = np.atleast_1d(np.asarray(self.tail.quantile(pt), dtype=float))
            else:
                out[hi] = [self._invert_cdf(float(pp)) for pp in ps[hi]]
        return float(out[0]) if arr.ndim == 0 else out.reshape(arr.shape)

    # --- sampling and moments ---
    def sample(self, size: int = 1) -> np.ndarray:
        if size <= 0:
            raise ValueError("size must be positive.")
        n_body = int(np.random.binomial(size, self.weight))
        n_tail = size - n_body
        parts = []
        if n_body > 0:
            u = np.random.uniform(0.0, self._Fb_u, size=n_body)
            parts.append(np.asarray(self.body.quantile(u), dtype=float).reshape(-1))
        if n_tail > 0:
            parts.append(np.asarray(self.tail.sample(n_tail), dtype=float).reshape(-1))
        out = np.concatenate(parts) if parts else np.empty(0)
        np.random.shuffle(out)
        return out

    def mean(self) -> float:
        if not hasattr(self.tail, "mean"):
            raise TypeError("tail must implement mean() for the spliced mean.")
        u, w, fb_u = self.threshold, self.weight, self._Fb_u
        lev = self.body.limited_expected_value(u)      # E[min(body, u)]
        body_integral = lev - u * (1.0 - fb_u)         # int_0^u t f_body(t) dt
        body_term = w * body_integral / fb_u
        tail_term = (1.0 - w) * float(self.tail.mean())
        return float(body_term + tail_term)

    def variance(self) -> float:
        if not (hasattr(self.tail, "mean") and hasattr(self.tail, "variance")):
            raise TypeError("tail must implement mean() and variance() for the spliced variance.")
        u, w, fb_u = self.threshold, self.weight, self._Fb_u
        body_second, _ = quad(lambda t: t * t * self.body.pdf(t), 0.0, u, limit=200)
        e_x2_body = w * body_second / fb_u
        tail_mean = float(self.tail.mean())
        e_x2_tail = (1.0 - w) * (float(self.tail.variance()) + tail_mean ** 2)
        m = self.mean()
        return float(e_x2_body + e_x2_tail - m * m)

    def __repr__(self) -> str:
        return (
            f"SplicedSeverity(body={self.body!r}, tail={self.tail!r}, "
            f"threshold={self.threshold}, weight={self.weight})"
        )

Weibull

Bases: SeverityModel

Weibull severity model.

Parameterization

X ~ Weibull(shape=k, scale=lam) Support: x > 0

Parameters

k : float Shape parameter, with k > 0. lam : float Scale parameter, with lam > 0.

Source code in lossmodels/severity/weibull.py 
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
class Weibull(SeverityModel):
    """
    Weibull severity model.

    Parameterization
    ----------------
    X ~ Weibull(shape=k, scale=lam)
    Support: x > 0

    Parameters
    ----------
    k : float
        Shape parameter, with k > 0.
    lam : float
        Scale parameter, with lam > 0.
    """

    def __init__(self, k: float, lam: float):
        if k <= 0:
            raise ValueError("k must be positive.")
        if lam <= 0:
            raise ValueError("lam must be positive.")

        self.k = k
        self.lam = lam

    def sample(self, size: int = 1) -> np.ndarray:
        if size <= 0:
            raise ValueError("size must be positive.")

        return self.lam * np.random.weibull(a=self.k, size=size)

    def mean(self) -> float:
        return float(self.lam * gamma_func(1 + 1 / self.k))

    def variance(self) -> float:
        m1 = gamma_func(1 + 1 / self.k)
        m2 = gamma_func(1 + 2 / self.k)
        return float(self.lam ** 2 * (m2 - m1 ** 2))

    def pdf(self, x):
        return eval_dist(lambda v: weibull_min.pdf(v, c=self.k, scale=self.lam), x)

    def cdf(self, x):
        return eval_dist(lambda v: weibull_min.cdf(v, c=self.k, scale=self.lam), x)

    def quantile(self, p):
        return eval_dist(lambda v: weibull_min.ppf(v, c=self.k, scale=self.lam), p)

    def __repr__(self) -> str:
        return f"Weibull(k={self.k}, lam={self.lam})"

Burr

Bases: SeverityModel

Burr (Type XII, Singh-Maddala) severity.

X ~ Burr(alpha, theta, gamma), support x > 0. F(x) = 1 - [1 / (1 + (x/theta)^gamma)]^alpha E[X^k] = theta^k Gamma(1 + k/gamma) Gamma(alpha - k/gamma) / Gamma(alpha), for -gamma < k < alpha * gamma.

Source code in lossmodels/severity/transformed_beta.py 
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
class Burr(SeverityModel):
    """Burr (Type XII, Singh-Maddala) severity.

    X ~ Burr(alpha, theta, gamma), support x > 0.
        F(x) = 1 - [1 / (1 + (x/theta)^gamma)]^alpha
        E[X^k] = theta^k Gamma(1 + k/gamma) Gamma(alpha - k/gamma) / Gamma(alpha),
                 for -gamma < k < alpha * gamma.
    """

    def __init__(self, alpha: float, theta: float, gamma: float):
        if alpha <= 0 or theta <= 0 or gamma <= 0:
            raise ValueError("alpha, theta, gamma must all be positive.")
        self.alpha = alpha
        self.theta = theta
        self.gamma = gamma
        self._d = burr12(c=gamma, d=alpha, scale=theta)

    def _moment(self, k: float) -> float:
        if not -self.gamma < k < self.alpha * self.gamma:
            raise ValueError(
                f"E[X^{k}] does not exist; need -gamma < k < alpha*gamma."
            )
        return (
            self.theta ** k
            * _G(1 + k / self.gamma)
            * _G(self.alpha - k / self.gamma)
            / _G(self.alpha)
        )

    def mean(self) -> float:
        if self.alpha * self.gamma <= 1:
            raise ValueError("Mean does not exist for alpha*gamma <= 1.")
        return self._moment(1)

    def variance(self) -> float:
        if self.alpha * self.gamma <= 2:
            raise ValueError("Variance does not exist for alpha*gamma <= 2.")
        return self._moment(2) - self._moment(1) ** 2

    def sample(self, size: int = 1) -> np.ndarray:
        if size <= 0:
            raise ValueError("size must be positive.")
        return self._d.rvs(size=size)

    def pdf(self, x):
        return eval_dist(lambda v: self._d.pdf(v), x)

    def cdf(self, x):
        return eval_dist(lambda v: self._d.cdf(v), x)

    def quantile(self, p):
        return eval_dist(lambda v: self._d.ppf(v), p)

    def __repr__(self) -> str:
        return f"Burr(alpha={self.alpha}, theta={self.theta}, gamma={self.gamma})"

InverseBurr

Bases: SeverityModel

Inverse Burr (Dagum) severity.

X ~ InverseBurr(tau, theta, gamma), support x > 0. F(x) = [ (x/theta)^gamma / (1 + (x/theta)^gamma) ]^tau E[X^k] = theta^k Gamma(tau + k/gamma) Gamma(1 - k/gamma) / Gamma(tau), for -tau*gamma < k < gamma.

Source code in lossmodels/severity/transformed_beta.py 
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
class InverseBurr(SeverityModel):
    """Inverse Burr (Dagum) severity.

    X ~ InverseBurr(tau, theta, gamma), support x > 0.
        F(x) = [ (x/theta)^gamma / (1 + (x/theta)^gamma) ]^tau
        E[X^k] = theta^k Gamma(tau + k/gamma) Gamma(1 - k/gamma) / Gamma(tau),
                 for -tau*gamma < k < gamma.
    """

    def __init__(self, tau: float, theta: float, gamma: float):
        if tau <= 0 or theta <= 0 or gamma <= 0:
            raise ValueError("tau, theta, gamma must all be positive.")
        self.tau = tau
        self.theta = theta
        self.gamma = gamma
        self._d = burr(c=gamma, d=tau, scale=theta)

    def _moment(self, k: float) -> float:
        if not -self.tau * self.gamma < k < self.gamma:
            raise ValueError(
                f"E[X^{k}] does not exist; need -tau*gamma < k < gamma."
            )
        return (
            self.theta ** k
            * _G(self.tau + k / self.gamma)
            * _G(1 - k / self.gamma)
            / _G(self.tau)
        )

    def mean(self) -> float:
        if self.gamma <= 1:
            raise ValueError("Mean does not exist for gamma <= 1.")
        return self._moment(1)

    def variance(self) -> float:
        if self.gamma <= 2:
            raise ValueError("Variance does not exist for gamma <= 2.")
        return self._moment(2) - self._moment(1) ** 2

    def sample(self, size: int = 1) -> np.ndarray:
        if size <= 0:
            raise ValueError("size must be positive.")
        return self._d.rvs(size=size)

    def pdf(self, x):
        return eval_dist(lambda v: self._d.pdf(v), x)

    def cdf(self, x):
        return eval_dist(lambda v: self._d.cdf(v), x)

    def quantile(self, p):
        return eval_dist(lambda v: self._d.ppf(v), p)

    def __repr__(self) -> str:
        return f"InverseBurr(tau={self.tau}, theta={self.theta}, gamma={self.gamma})"

GeneralizedPareto

Bases: SeverityModel

Generalized Pareto severity (Klugman transformed-beta form; NOT the EVT GPD).

X ~ GeneralizedPareto(alpha, theta, tau), support x > 0. F(x) = Beta(tau, alpha; x/(x+theta)) (regularized incomplete beta) E[X^k] = theta^k Gamma(tau + k) Gamma(alpha - k) / (Gamma(alpha) Gamma(tau)), for -tau < k < alpha.

The extreme-value GPD (peaks-over-threshold tail) lives in extremeloss; this is the three-parameter loss-severity distribution from Loss Models.

Source code in lossmodels/severity/transformed_beta.py 
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
class GeneralizedPareto(SeverityModel):
    """Generalized Pareto severity (Klugman transformed-beta form; NOT the EVT GPD).

    X ~ GeneralizedPareto(alpha, theta, tau), support x > 0.
        F(x) = Beta(tau, alpha; x/(x+theta))   (regularized incomplete beta)
        E[X^k] = theta^k Gamma(tau + k) Gamma(alpha - k) / (Gamma(alpha) Gamma(tau)),
                 for -tau < k < alpha.

    The extreme-value GPD (peaks-over-threshold tail) lives in ``extremeloss``;
    this is the three-parameter loss-severity distribution from Loss Models.
    """

    def __init__(self, alpha: float, theta: float, tau: float):
        if alpha <= 0 or theta <= 0 or tau <= 0:
            raise ValueError("alpha, theta, tau must all be positive.")
        self.alpha = alpha
        self.theta = theta
        self.tau = tau
        self._d = betaprime(a=tau, b=alpha, scale=theta)

    def _moment(self, k: float) -> float:
        if not -self.tau < k < self.alpha:
            raise ValueError(f"E[X^{k}] does not exist; need -tau < k < alpha.")
        return (
            self.theta ** k
            * _G(self.tau + k)
            * _G(self.alpha - k)
            / (_G(self.alpha) * _G(self.tau))
        )

    def mean(self) -> float:
        if self.alpha <= 1:
            raise ValueError("Mean does not exist for alpha <= 1.")
        return self._moment(1)

    def variance(self) -> float:
        if self.alpha <= 2:
            raise ValueError("Variance does not exist for alpha <= 2.")
        return self._moment(2) - self._moment(1) ** 2

    def sample(self, size: int = 1) -> np.ndarray:
        if size <= 0:
            raise ValueError("size must be positive.")
        return self._d.rvs(size=size)

    def pdf(self, x):
        return eval_dist(lambda v: self._d.pdf(v), x)

    def cdf(self, x):
        return eval_dist(lambda v: self._d.cdf(v), x)

    def quantile(self, p):
        return eval_dist(lambda v: self._d.ppf(v), p)

    def __repr__(self) -> str:
        return (
            f"GeneralizedPareto(alpha={self.alpha}, theta={self.theta}, "
            f"tau={self.tau})"
        )

ParetoII

Bases: SeverityModel

Pareto (Type II, Lomax) severity -- the FAM/ASTAM two-parameter "Pareto".

X ~ Pareto(alpha, theta), support x > 0. f(x) = alpha theta^alpha / (x + theta)^(alpha+1) F(x) = 1 - [theta / (x + theta)]^alpha E[X^k] = theta^k k! / [(alpha-1)...(alpha-k)] (positive integer k < alpha) E[X] = theta / (alpha - 1), alpha > 1.

Note: this is distinct from :class:lossmodels.Pareto, which is the Pareto Type I (single-parameter) distribution with support x >= theta. The FAM table lists this Type II / Lomax form under the name "Pareto"; the Type I appears under "Single-Parameter Pareto".

Source code in lossmodels/severity/transformed_beta.py 
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
class ParetoII(SeverityModel):
    """Pareto (Type II, Lomax) severity -- the FAM/ASTAM two-parameter "Pareto".

    X ~ Pareto(alpha, theta), support x > 0.
        f(x) = alpha theta^alpha / (x + theta)^(alpha+1)
        F(x) = 1 - [theta / (x + theta)]^alpha
        E[X^k] = theta^k k! / [(alpha-1)...(alpha-k)]  (positive integer k < alpha)
        E[X] = theta / (alpha - 1),  alpha > 1.

    Note: this is distinct from :class:`lossmodels.Pareto`, which is the Pareto
    Type I (single-parameter) distribution with support x >= theta. The FAM table
    lists this Type II / Lomax form under the name "Pareto"; the Type I appears
    under "Single-Parameter Pareto".
    """

    def __init__(self, alpha: float, theta: float):
        if alpha <= 0 or theta <= 0:
            raise ValueError("alpha and theta must be positive.")
        self.alpha = alpha
        self.theta = theta
        self._d = lomax(c=alpha, scale=theta)

    def _moment(self, k: float) -> float:
        if not -1 < k < self.alpha:
            raise ValueError(f"E[X^{k}] does not exist; need -1 < k < alpha.")
        return self.theta ** k * _G(k + 1) * _G(self.alpha - k) / _G(self.alpha)

    def mean(self) -> float:
        if self.alpha <= 1:
            raise ValueError("Mean does not exist for alpha <= 1.")
        return self.theta / (self.alpha - 1)

    def variance(self) -> float:
        if self.alpha <= 2:
            raise ValueError("Variance does not exist for alpha <= 2.")
        return self._moment(2) - self._moment(1) ** 2

    def sample(self, size: int = 1) -> np.ndarray:
        if size <= 0:
            raise ValueError("size must be positive.")
        return self._d.rvs(size=size)

    def pdf(self, x):
        return eval_dist(lambda v: self._d.pdf(v), x)

    def cdf(self, x):
        return eval_dist(lambda v: self._d.cdf(v), x)

    def quantile(self, p):
        return eval_dist(lambda v: self._d.ppf(v), p)

    def __repr__(self) -> str:
        return f"ParetoII(alpha={self.alpha}, theta={self.theta})"

InversePareto

Bases: SeverityModel

Inverse Pareto severity.

X ~ InversePareto(tau, theta), support x > 0. f(x) = tau theta x^(tau-1) / (x + theta)^(tau+1) F(x) = [x / (x + theta)]^tau E[X^k] = theta^k Gamma(tau + k) Gamma(1 - k) / Gamma(tau), -tau < k < 1.

The mean (k = 1) lies on the boundary of the moment range and does not exist, so :meth:mean and :meth:variance always raise.

Source code in lossmodels/severity/transformed_beta.py 
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
class InversePareto(SeverityModel):
    """Inverse Pareto severity.

    X ~ InversePareto(tau, theta), support x > 0.
        f(x) = tau theta x^(tau-1) / (x + theta)^(tau+1)
        F(x) = [x / (x + theta)]^tau
        E[X^k] = theta^k Gamma(tau + k) Gamma(1 - k) / Gamma(tau),  -tau < k < 1.

    The mean (k = 1) lies on the boundary of the moment range and does not exist,
    so :meth:`mean` and :meth:`variance` always raise.
    """

    def __init__(self, tau: float, theta: float):
        if tau <= 0 or theta <= 0:
            raise ValueError("tau and theta must be positive.")
        self.tau = tau
        self.theta = theta
        self._d = betaprime(a=tau, b=1, scale=theta)

    def _moment(self, k: float) -> float:
        if not -self.tau < k < 1:
            raise ValueError(f"E[X^{k}] does not exist; need -tau < k < 1.")
        return self.theta ** k * _G(self.tau + k) * _G(1 - k) / _G(self.tau)

    def mean(self) -> float:
        raise ValueError("Mean does not exist for the inverse Pareto distribution.")

    def variance(self) -> float:
        raise ValueError(
            "Variance does not exist for the inverse Pareto distribution."
        )

    def sample(self, size: int = 1) -> np.ndarray:
        if size <= 0:
            raise ValueError("size must be positive.")
        return self._d.rvs(size=size)

    def pdf(self, x):
        return eval_dist(lambda v: self._d.pdf(v), x)

    def cdf(self, x):
        return eval_dist(lambda v: self._d.cdf(v), x)

    def quantile(self, p):
        return eval_dist(lambda v: self._d.ppf(v), p)

    def __repr__(self) -> str:
        return f"InversePareto(tau={self.tau}, theta={self.theta})"

Loglogistic

Bases: SeverityModel

Loglogistic (Fisk) severity.

X ~ Loglogistic(gamma, theta), support x > 0. F(x) = (x/theta)^gamma / (1 + (x/theta)^gamma) E[X^k] = theta^k Gamma(1 + k/gamma) Gamma(1 - k/gamma), -gamma < k < gamma.

Source code in lossmodels/severity/transformed_beta.py 
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
class Loglogistic(SeverityModel):
    """Loglogistic (Fisk) severity.

    X ~ Loglogistic(gamma, theta), support x > 0.
        F(x) = (x/theta)^gamma / (1 + (x/theta)^gamma)
        E[X^k] = theta^k Gamma(1 + k/gamma) Gamma(1 - k/gamma),  -gamma < k < gamma.
    """

    def __init__(self, gamma: float, theta: float):
        if gamma <= 0 or theta <= 0:
            raise ValueError("gamma and theta must be positive.")
        self.gamma = gamma
        self.theta = theta
        self._d = fisk(c=gamma, scale=theta)

    def _moment(self, k: float) -> float:
        if not -self.gamma < k < self.gamma:
            raise ValueError(f"E[X^{k}] does not exist; need -gamma < k < gamma.")
        return self.theta ** k * _G(1 + k / self.gamma) * _G(1 - k / self.gamma)

    def mean(self) -> float:
        if self.gamma <= 1:
            raise ValueError("Mean does not exist for gamma <= 1.")
        return self._moment(1)

    def variance(self) -> float:
        if self.gamma <= 2:
            raise ValueError("Variance does not exist for gamma <= 2.")
        return self._moment(2) - self._moment(1) ** 2

    def sample(self, size: int = 1) -> np.ndarray:
        if size <= 0:
            raise ValueError("size must be positive.")
        return self._d.rvs(size=size)

    def pdf(self, x):
        return eval_dist(lambda v: self._d.pdf(v), x)

    def cdf(self, x):
        return eval_dist(lambda v: self._d.cdf(v), x)

    def quantile(self, p):
        return eval_dist(lambda v: self._d.ppf(v), p)

    def __repr__(self) -> str:
        return f"Loglogistic(gamma={self.gamma}, theta={self.theta})"

Paralogistic

Bases: SeverityModel

Paralogistic severity (a Burr distribution with gamma = alpha).

X ~ Paralogistic(alpha, theta), support x > 0. E[X^k] = theta^k Gamma(1 + k/alpha) Gamma(alpha - k/alpha) / Gamma(alpha), for -alpha < k < alpha^2.

Source code in lossmodels/severity/transformed_beta.py 
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
class Paralogistic(SeverityModel):
    """Paralogistic severity (a Burr distribution with gamma = alpha).

    X ~ Paralogistic(alpha, theta), support x > 0.
        E[X^k] = theta^k Gamma(1 + k/alpha) Gamma(alpha - k/alpha) / Gamma(alpha),
                 for -alpha < k < alpha^2.
    """

    def __init__(self, alpha: float, theta: float):
        if alpha <= 0 or theta <= 0:
            raise ValueError("alpha and theta must be positive.")
        self.alpha = alpha
        self.theta = theta
        self._d = burr12(c=alpha, d=alpha, scale=theta)

    def _moment(self, k: float) -> float:
        if not -self.alpha < k < self.alpha ** 2:
            raise ValueError(
                f"E[X^{k}] does not exist; need -alpha < k < alpha^2."
            )
        return (
            self.theta ** k
            * _G(1 + k / self.alpha)
            * _G(self.alpha - k / self.alpha)
            / _G(self.alpha)
        )

    def mean(self) -> float:
        if self.alpha ** 2 <= 1:
            raise ValueError("Mean does not exist for alpha^2 <= 1.")
        return self._moment(1)

    def variance(self) -> float:
        if self.alpha ** 2 <= 2:
            raise ValueError("Variance does not exist for alpha^2 <= 2.")
        return self._moment(2) - self._moment(1) ** 2

    def sample(self, size: int = 1) -> np.ndarray:
        if size <= 0:
            raise ValueError("size must be positive.")
        return self._d.rvs(size=size)

    def pdf(self, x):
        return eval_dist(lambda v: self._d.pdf(v), x)

    def cdf(self, x):
        return eval_dist(lambda v: self._d.cdf(v), x)

    def quantile(self, p):
        return eval_dist(lambda v: self._d.ppf(v), p)

    def __repr__(self) -> str:
        return f"Paralogistic(alpha={self.alpha}, theta={self.theta})"

InverseParalogistic

Bases: SeverityModel

Inverse paralogistic severity (an inverse Burr with gamma = tau).

X ~ InverseParalogistic(tau, theta), support x > 0. E[X^k] = theta^k Gamma(tau + k/tau) Gamma(1 - k/tau) / Gamma(tau), for -tau^2 < k < tau.

Source code in lossmodels/severity/transformed_beta.py 
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
class InverseParalogistic(SeverityModel):
    """Inverse paralogistic severity (an inverse Burr with gamma = tau).

    X ~ InverseParalogistic(tau, theta), support x > 0.
        E[X^k] = theta^k Gamma(tau + k/tau) Gamma(1 - k/tau) / Gamma(tau),
                 for -tau^2 < k < tau.
    """

    def __init__(self, tau: float, theta: float):
        if tau <= 0 or theta <= 0:
            raise ValueError("tau and theta must be positive.")
        self.tau = tau
        self.theta = theta
        self._d = burr(c=tau, d=tau, scale=theta)

    def _moment(self, k: float) -> float:
        if not -self.tau ** 2 < k < self.tau:
            raise ValueError(f"E[X^{k}] does not exist; need -tau^2 < k < tau.")
        return (
            self.theta ** k
            * _G(self.tau + k / self.tau)
            * _G(1 - k / self.tau)
            / _G(self.tau)
        )

    def mean(self) -> float:
        if self.tau <= 1:
            raise ValueError("Mean does not exist for tau <= 1.")
        return self._moment(1)

    def variance(self) -> float:
        if self.tau <= 2:
            raise ValueError("Variance does not exist for tau <= 2.")
        return self._moment(2) - self._moment(1) ** 2

    def sample(self, size: int = 1) -> np.ndarray:
        if size <= 0:
            raise ValueError("size must be positive.")
        return self._d.rvs(size=size)

    def pdf(self, x):
        return eval_dist(lambda v: self._d.pdf(v), x)

    def cdf(self, x):
        return eval_dist(lambda v: self._d.cdf(v), x)

    def quantile(self, p):
        return eval_dist(lambda v: self._d.ppf(v), p)

    def __repr__(self) -> str:
        return f"InverseParalogistic(tau={self.tau}, theta={self.theta})"

InverseGamma

Bases: SeverityModel

Inverse gamma (Vinci) severity.

X ~ InverseGamma(alpha, theta), support x > 0. F(x) = 1 - Gamma(alpha; theta/x) E[X^k] = theta^k Gamma(alpha - k) / Gamma(alpha), k < alpha.

Source code in lossmodels/severity/transformed_gamma.py 
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
class InverseGamma(SeverityModel):
    """Inverse gamma (Vinci) severity.

    X ~ InverseGamma(alpha, theta), support x > 0.
        F(x) = 1 - Gamma(alpha; theta/x)
        E[X^k] = theta^k Gamma(alpha - k) / Gamma(alpha),  k < alpha.
    """

    def __init__(self, alpha: float, theta: float):
        if alpha <= 0 or theta <= 0:
            raise ValueError("alpha and theta must be positive.")
        self.alpha = alpha
        self.theta = theta
        self._d = invgamma(a=alpha, scale=theta)

    def _moment(self, k: float) -> float:
        if k >= self.alpha:
            raise ValueError(f"E[X^{k}] does not exist; need k < alpha.")
        return self.theta ** k * _G(self.alpha - k) / _G(self.alpha)

    def mean(self) -> float:
        if self.alpha <= 1:
            raise ValueError("Mean does not exist for alpha <= 1.")
        return self._moment(1)

    def variance(self) -> float:
        if self.alpha <= 2:
            raise ValueError("Variance does not exist for alpha <= 2.")
        return self._moment(2) - self._moment(1) ** 2

    def sample(self, size: int = 1) -> np.ndarray:
        if size <= 0:
            raise ValueError("size must be positive.")
        return self._d.rvs(size=size)

    def pdf(self, x):
        return eval_dist(lambda v: self._d.pdf(v), x)

    def cdf(self, x):
        return eval_dist(lambda v: self._d.cdf(v), x)

    def quantile(self, p):
        return eval_dist(lambda v: self._d.ppf(v), p)

    def __repr__(self) -> str:
        return f"InverseGamma(alpha={self.alpha}, theta={self.theta})"

InverseWeibull

Bases: SeverityModel

Inverse Weibull (log-Gompertz) severity.

X ~ InverseWeibull(theta, tau), support x > 0. F(x) = exp[-(theta/x)^tau] E[X^k] = theta^k Gamma(1 - k/tau), k < tau.

Source code in lossmodels/severity/transformed_gamma.py 
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
class InverseWeibull(SeverityModel):
    """Inverse Weibull (log-Gompertz) severity.

    X ~ InverseWeibull(theta, tau), support x > 0.
        F(x) = exp[-(theta/x)^tau]
        E[X^k] = theta^k Gamma(1 - k/tau),  k < tau.
    """

    def __init__(self, theta: float, tau: float):
        if theta <= 0 or tau <= 0:
            raise ValueError("theta and tau must be positive.")
        self.theta = theta
        self.tau = tau
        self._d = invweibull(c=tau, scale=theta)

    def _moment(self, k: float) -> float:
        if k >= self.tau:
            raise ValueError(f"E[X^{k}] does not exist; need k < tau.")
        return self.theta ** k * _G(1 - k / self.tau)

    def mean(self) -> float:
        if self.tau <= 1:
            raise ValueError("Mean does not exist for tau <= 1.")
        return self._moment(1)

    def variance(self) -> float:
        if self.tau <= 2:
            raise ValueError("Variance does not exist for tau <= 2.")
        return self._moment(2) - self._moment(1) ** 2

    def sample(self, size: int = 1) -> np.ndarray:
        if size <= 0:
            raise ValueError("size must be positive.")
        return self._d.rvs(size=size)

    def pdf(self, x):
        return eval_dist(lambda v: self._d.pdf(v), x)

    def cdf(self, x):
        return eval_dist(lambda v: self._d.cdf(v), x)

    def quantile(self, p):
        return eval_dist(lambda v: self._d.ppf(v), p)

    def __repr__(self) -> str:
        return f"InverseWeibull(theta={self.theta}, tau={self.tau})"

InverseExponential

Bases: SeverityModel

Inverse exponential severity.

X ~ InverseExponential(theta), support x > 0. F(x) = exp(-theta/x) E[X^k] = theta^k Gamma(1 - k), k < 1.

The mean (k = 1) does not exist, so :meth:mean and :meth:variance raise.

Source code in lossmodels/severity/transformed_gamma.py 
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
class InverseExponential(SeverityModel):
    """Inverse exponential severity.

    X ~ InverseExponential(theta), support x > 0.
        F(x) = exp(-theta/x)
        E[X^k] = theta^k Gamma(1 - k),  k < 1.

    The mean (k = 1) does not exist, so :meth:`mean` and :meth:`variance` raise.
    """

    def __init__(self, theta: float):
        if theta <= 0:
            raise ValueError("theta must be positive.")
        self.theta = theta
        self._d = invweibull(c=1, scale=theta)

    def _moment(self, k: float) -> float:
        if k >= 1:
            raise ValueError(f"E[X^{k}] does not exist; need k < 1.")
        return self.theta ** k * _G(1 - k)

    def mean(self) -> float:
        raise ValueError(
            "Mean does not exist for the inverse exponential distribution."
        )

    def variance(self) -> float:
        raise ValueError(
            "Variance does not exist for the inverse exponential distribution."
        )

    def sample(self, size: int = 1) -> np.ndarray:
        if size <= 0:
            raise ValueError("size must be positive.")
        return self._d.rvs(size=size)

    def pdf(self, x):
        return eval_dist(lambda v: self._d.pdf(v), x)

    def cdf(self, x):
        return eval_dist(lambda v: self._d.cdf(v), x)

    def quantile(self, p):
        return eval_dist(lambda v: self._d.ppf(v), p)

    def __repr__(self) -> str:
        return f"InverseExponential(theta={self.theta})"

InverseGaussian

Bases: SeverityModel

Inverse Gaussian (Wald) severity, Klugman (mu, theta) parameterization.

X ~ InverseGaussian(mu, theta), support x > 0. E[X] = mu, Var[X] = mu^3 / theta.

Mapped to scipy.stats.invgauss(mu=mu/theta, scale=theta).

Source code in lossmodels/severity/other_severity.py 
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
class InverseGaussian(SeverityModel):
    """Inverse Gaussian (Wald) severity, Klugman ``(mu, theta)`` parameterization.

    X ~ InverseGaussian(mu, theta), support x > 0.
        E[X] = mu,  Var[X] = mu^3 / theta.

    Mapped to ``scipy.stats.invgauss(mu=mu/theta, scale=theta)``.
    """

    def __init__(self, mu: float, theta: float):
        if mu <= 0 or theta <= 0:
            raise ValueError("mu and theta must be positive.")
        self.mu = mu
        self.theta = theta
        self._d = invgauss(mu=mu / theta, scale=theta)

    def mean(self) -> float:
        return self.mu

    def variance(self) -> float:
        return self.mu ** 3 / self.theta

    def sample(self, size: int = 1) -> np.ndarray:
        if size <= 0:
            raise ValueError("size must be positive.")
        return self._d.rvs(size=size)

    def pdf(self, x):
        return eval_dist(lambda v: self._d.pdf(v), x)

    def cdf(self, x):
        return eval_dist(lambda v: self._d.cdf(v), x)

    def quantile(self, p):
        return eval_dist(lambda v: self._d.ppf(v), p)

    def __repr__(self) -> str:
        return f"InverseGaussian(mu={self.mu}, theta={self.theta})"

LogT

Bases: SeverityModel

Log-t severity.

If Y has a Student-t distribution with r degrees of freedom, then X = exp(sigma * Y + mu) has the log-t distribution (support x > 0). It is heavier-tailed than the lognormal; no positive moments exist, so :meth:mean and :meth:variance raise. mu may be negative.

F(x) = F_r((ln x - mu) / sigma),   F_r the Student-t cdf.
Source code in lossmodels/severity/other_severity.py 
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
class LogT(SeverityModel):
    """Log-t severity.

    If Y has a Student-t distribution with r degrees of freedom, then
    X = exp(sigma * Y + mu) has the log-t distribution (support x > 0). It is
    heavier-tailed than the lognormal; **no positive moments exist**, so
    :meth:`mean` and :meth:`variance` raise. ``mu`` may be negative.

        F(x) = F_r((ln x - mu) / sigma),   F_r the Student-t cdf.
    """

    def __init__(self, r: float, mu: float, sigma: float):
        if r <= 0:
            raise ValueError("r (degrees of freedom) must be positive.")
        if sigma <= 0:
            raise ValueError("sigma must be positive.")
        self.r = r
        self.mu = mu
        self.sigma = sigma
        self._t = student_t(df=r)

    def mean(self) -> float:
        raise ValueError("Positive moments do not exist for the log-t distribution.")

    def variance(self) -> float:
        raise ValueError("Positive moments do not exist for the log-t distribution.")

    def sample(self, size: int = 1) -> np.ndarray:
        if size <= 0:
            raise ValueError("size must be positive.")
        return np.exp(self.sigma * self._t.rvs(size=size) + self.mu)

    def pdf(self, x):
        def f(v):
            v = np.asarray(v, dtype=float)
            pos = v > 0
            safe = np.where(pos, v, 1.0)
            z = (np.log(safe) - self.mu) / self.sigma
            return np.where(pos, self._t.pdf(z) / (safe * self.sigma), 0.0)

        return eval_dist(f, x)

    def cdf(self, x):
        def f(v):
            v = np.asarray(v, dtype=float)
            pos = v > 0
            safe = np.where(pos, v, 1.0)
            z = (np.log(safe) - self.mu) / self.sigma
            return np.where(pos, self._t.cdf(z), 0.0)

        return eval_dist(f, x)

    def quantile(self, p):
        return eval_dist(
            lambda v: np.exp(self.sigma * self._t.ppf(v) + self.mu), p
        )

    def __repr__(self) -> str:
        return f"LogT(r={self.r}, mu={self.mu}, sigma={self.sigma})"

SingleParameterPareto

Bases: SeverityModel

Single-Parameter Pareto severity (Pareto Type I; theta is a fixed lower bound).

X ~ SingleParameterPareto(alpha, theta), support x > theta. f(x) = alpha theta^alpha / x^(alpha+1), x > theta F(x) = 1 - (theta/x)^alpha, x > theta E[X^k] = alpha theta^k / (alpha - k), k < alpha.

This is the same distribution as :class:lossmodels.Pareto; it is provided under the FAM/ASTAM name. Per Loss Models, only alpha is a true parameter -- theta is the known lower truncation point set in advance.

Source code in lossmodels/severity/other_severity.py 
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
class SingleParameterPareto(SeverityModel):
    """Single-Parameter Pareto severity (Pareto Type I; theta is a fixed lower bound).

    X ~ SingleParameterPareto(alpha, theta), support x > theta.
        f(x) = alpha theta^alpha / x^(alpha+1),  x > theta
        F(x) = 1 - (theta/x)^alpha,  x > theta
        E[X^k] = alpha theta^k / (alpha - k),  k < alpha.

    This is the same distribution as :class:`lossmodels.Pareto`; it is provided
    under the FAM/ASTAM name. Per Loss Models, only ``alpha`` is a true parameter
    -- ``theta`` is the known lower truncation point set in advance.
    """

    def __init__(self, alpha: float, theta: float):
        if alpha <= 0 or theta <= 0:
            raise ValueError("alpha and theta must be positive.")
        self.alpha = alpha
        self.theta = theta
        self._d = pareto(b=alpha, scale=theta)

    def _moment(self, k: float) -> float:
        if k >= self.alpha:
            raise ValueError(f"E[X^{k}] does not exist; need k < alpha.")
        return self.alpha * self.theta ** k / (self.alpha - k)

    def mean(self) -> float:
        if self.alpha <= 1:
            raise ValueError("Mean does not exist for alpha <= 1.")
        return self._moment(1)

    def variance(self) -> float:
        if self.alpha <= 2:
            raise ValueError("Variance does not exist for alpha <= 2.")
        return self._moment(2) - self._moment(1) ** 2

    def sample(self, size: int = 1) -> np.ndarray:
        if size <= 0:
            raise ValueError("size must be positive.")
        return self._d.rvs(size=size)

    def pdf(self, x):
        return eval_dist(lambda v: self._d.pdf(v), x)

    def cdf(self, x):
        return eval_dist(lambda v: self._d.cdf(v), x)

    def quantile(self, p):
        return eval_dist(lambda v: self._d.ppf(v), p)

    def __repr__(self) -> str:
        return f"SingleParameterPareto(alpha={self.alpha}, theta={self.theta})"

Beta

Bases: SeverityModel

Beta severity on (0, theta).

X ~ Beta(a, b, theta), support 0 < x < theta. E[X^k] = theta^k Gamma(a+b) Gamma(a+k) / (Gamma(a) Gamma(a+b+k)), k > -a. E[X] = theta a / (a + b).

Source code in lossmodels/severity/finite_support.py 
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
class Beta(SeverityModel):
    """Beta severity on (0, theta).

    X ~ Beta(a, b, theta), support 0 < x < theta.
        E[X^k] = theta^k Gamma(a+b) Gamma(a+k) / (Gamma(a) Gamma(a+b+k)),  k > -a.
        E[X] = theta a / (a + b).
    """

    def __init__(self, a: float, b: float, theta: float):
        if a <= 0 or b <= 0 or theta <= 0:
            raise ValueError("a, b, theta must all be positive.")
        self.a = a
        self.b = b
        self.theta = theta
        self._d = beta_dist(a, b, scale=theta)

    def _moment(self, k: float) -> float:
        return (
            self.theta ** k
            * _G(self.a + self.b)
            * _G(self.a + k)
            / (_G(self.a) * _G(self.a + self.b + k))
        )

    def mean(self) -> float:
        return self.theta * self.a / (self.a + self.b)

    def variance(self) -> float:
        return self._moment(2) - self._moment(1) ** 2

    def sample(self, size: int = 1) -> np.ndarray:
        if size <= 0:
            raise ValueError("size must be positive.")
        return self._d.rvs(size=size)

    def pdf(self, x):
        return eval_dist(lambda v: self._d.pdf(v), x)

    def cdf(self, x):
        return eval_dist(lambda v: self._d.cdf(v), x)

    def quantile(self, p):
        return eval_dist(lambda v: self._d.ppf(v), p)

    def __repr__(self) -> str:
        return f"Beta(a={self.a}, b={self.b}, theta={self.theta})"

GeneralizedBeta

Bases: SeverityModel

Generalized beta severity on (0, theta).

X ~ GeneralizedBeta(a, b, theta, tau), support 0 < x < theta, defined by U = (X/theta)^tau ~ Beta(a, b). Equivalently X = theta * U^(1/tau). F(x) = Beta(a, b; (x/theta)^tau) E[X^k] = theta^k Gamma(a+b) Gamma(a + k/tau) / (Gamma(a) Gamma(a + b + k/tau)), k > -a*tau.

Source code in lossmodels/severity/finite_support.py 
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
class GeneralizedBeta(SeverityModel):
    """Generalized beta severity on (0, theta).

    X ~ GeneralizedBeta(a, b, theta, tau), support 0 < x < theta, defined by
    U = (X/theta)^tau ~ Beta(a, b). Equivalently X = theta * U^(1/tau).
        F(x) = Beta(a, b; (x/theta)^tau)
        E[X^k] = theta^k Gamma(a+b) Gamma(a + k/tau) / (Gamma(a) Gamma(a + b + k/tau)),
                 k > -a*tau.
    """

    def __init__(self, a: float, b: float, theta: float, tau: float):
        if a <= 0 or b <= 0 or theta <= 0 or tau <= 0:
            raise ValueError("a, b, theta, tau must all be positive.")
        self.a = a
        self.b = b
        self.theta = theta
        self.tau = tau
        self._b = beta_dist(a, b)  # standard Beta(a, b) on (0, 1)

    def _moment(self, k: float) -> float:
        return (
            self.theta ** k
            * _G(self.a + self.b)
            * _G(self.a + k / self.tau)
            / (_G(self.a) * _G(self.a + self.b + k / self.tau))
        )

    def mean(self) -> float:
        return self._moment(1)

    def variance(self) -> float:
        return self._moment(2) - self._moment(1) ** 2

    def sample(self, size: int = 1) -> np.ndarray:
        if size <= 0:
            raise ValueError("size must be positive.")
        return self.theta * self._b.rvs(size=size) ** (1.0 / self.tau)

    def pdf(self, x):
        def f(v):
            v = np.asarray(v, dtype=float)
            pos = v > 0
            safe = np.where(pos, v, 1.0)
            u = (safe / self.theta) ** self.tau
            dens = (
                self._b.pdf(u)
                * self.tau
                * safe ** (self.tau - 1)
                / self.theta ** self.tau
            )
            return np.where(pos, dens, 0.0)

        return eval_dist(f, x)

    def cdf(self, x):
        def f(v):
            v = np.asarray(v, dtype=float)
            pos = v > 0
            safe = np.where(pos, v, 1.0)
            u = np.clip((safe / self.theta) ** self.tau, 0.0, 1.0)
            return np.where(v <= 0, 0.0, self._b.cdf(u))

        return eval_dist(f, x)

    def quantile(self, p):
        return eval_dist(
            lambda v: self.theta * self._b.ppf(v) ** (1.0 / self.tau), p
        )

    def __repr__(self) -> str:
        return (
            f"GeneralizedBeta(a={self.a}, b={self.b}, theta={self.theta}, "
            f"tau={self.tau})"
        )

FrequencyModel

Bases: ABC

Base class for all frequency distributions.

All frequency models must implement: - sample - mean - variance

Source code in lossmodels/frequency/base.py 
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
class FrequencyModel(ABC):
    """
    Base class for all frequency distributions.

    All frequency models must implement:
    - sample
    - mean
    - variance
    """

    @abstractmethod
    def sample(self, size: int = 1) -> np.ndarray:
        """
        Generate random samples of claim counts.

        Parameters
        ----------
        size : int
            Number of samples

        Returns
        -------
        np.ndarray
            Array of claim counts
        """
        pass

    @abstractmethod
    def mean(self) -> float:
        """Expected number of claims"""
        pass

    @abstractmethod
    def variance(self) -> float:
        """Variance of number of claims"""
        pass

    def std(self) -> float:
        """Standard deviation of claim count"""
        return np.sqrt(self.variance())

    def __repr__(self):
        return f"{self.__class__.__name__}()"
sample abstractmethod 
sample(size: int = 1) -> np.ndarray

Generate random samples of claim counts.

Parameters

size : int Number of samples

Returns

np.ndarray Array of claim counts

Source code in lossmodels/frequency/base.py 
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
@abstractmethod
def sample(self, size: int = 1) -> np.ndarray:
    """
    Generate random samples of claim counts.

    Parameters
    ----------
    size : int
        Number of samples

    Returns
    -------
    np.ndarray
        Array of claim counts
    """
    pass
mean abstractmethod 
mean() -> float

Expected number of claims

Source code in lossmodels/frequency/base.py 
32
33
34
35
@abstractmethod
def mean(self) -> float:
    """Expected number of claims"""
    pass
variance abstractmethod 
variance() -> float

Variance of number of claims

Source code in lossmodels/frequency/base.py 
37
38
39
40
@abstractmethod
def variance(self) -> float:
    """Variance of number of claims"""
    pass
std 
std() -> float

Standard deviation of claim count

Source code in lossmodels/frequency/base.py 
42
43
44
def std(self) -> float:
    """Standard deviation of claim count"""
    return np.sqrt(self.variance())

Binomial

Bases: FrequencyModel

Binomial frequency model.

Parameters

n : int Number of trials p : float Probability of success

Source code in lossmodels/frequency/binomial.py 
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
class Binomial(FrequencyModel):
    """
    Binomial frequency model.

    Parameters
    ----------
    n : int
        Number of trials
    p : float
        Probability of success
    """

    def __init__(self, n: int, p: float):
        if n <= 0:
            raise ValueError("n must be positive")
        if not (0 <= p <= 1):
            raise ValueError("p must be between 0 and 1")

        self.n = n
        self.p = p

    def sample(self, size: int = 1) -> np.ndarray:
        return np.random.binomial(self.n, self.p, size=size)

    def mean(self) -> float:
        return self.n * self.p

    def variance(self) -> float:
        return self.n * self.p * (1 - self.p)

    def pmf(self, k):
        """Probability mass function P(N = k)."""
        return eval_dist(lambda v: binom.pmf(v, self.n, self.p), k)

    def cdf(self, k):
        """Cumulative distribution function P(N <= k)."""
        return eval_dist(lambda v: binom.cdf(v, self.n, self.p), k)

    def __repr__(self):
        return f"Binomial(n={self.n}, p={self.p})"
pmf 
pmf(k)

Probability mass function P(N = k).

Source code in lossmodels/frequency/binomial.py 
38
39
40
def pmf(self, k):
    """Probability mass function P(N = k)."""
    return eval_dist(lambda v: binom.pmf(v, self.n, self.p), k)
cdf 
cdf(k)

Cumulative distribution function P(N <= k).

Source code in lossmodels/frequency/binomial.py 
42
43
44
def cdf(self, k):
    """Cumulative distribution function P(N <= k)."""
    return eval_dist(lambda v: binom.cdf(v, self.n, self.p), k)

Geometric

Bases: FrequencyModel

Geometric frequency model.

Support starting at 0: {0, 1, 2, 3, ...}

Parameters

p : float Probability of success

Notes

NumPy and SciPy define the geometric distribution on {1, 2, 3, ...} as the number of trials until first success. This implementation shifts that convention by 1 so the support starts at 0, which is more natural for claim counts.

Source code in lossmodels/frequency/geometric.py 
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
class Geometric(FrequencyModel):
    """
    Geometric frequency model.

    Support starting at 0: {0, 1, 2, 3, ...}

    Parameters
    ----------
    p : float
        Probability of success

    Notes
    -----
    NumPy and SciPy define the geometric distribution on {1, 2, 3, ...}
    as the number of trials until first success. This implementation shifts
    that convention by 1 so the support starts at 0, which is more natural
    for claim counts.
    """

    def __init__(self, p: float):
        if not (0 < p <= 1):
            raise ValueError("p must be in (0, 1].")

        self.p = p

    def sample(self, size: int = 1) -> np.ndarray:
        return np.random.geometric(self.p, size=size) - 1

    def mean(self) -> float:
        return (1 - self.p) / self.p

    def variance(self) -> float:
        return (1 - self.p) / (self.p ** 2)

    def pmf(self, k):
        """Probability mass function P(N = k)."""
        return eval_dist(lambda v: geom.pmf(v + 1, self.p), k)

    def cdf(self, k):
        """Cumulative distribution function P(N <= k)."""
        return eval_dist(lambda v: geom.cdf(v + 1, self.p), k)

    def __repr__(self):
        return f"Geometric(p={self.p})"
pmf 
pmf(k)

Probability mass function P(N = k).

Source code in lossmodels/frequency/geometric.py 
42
43
44
def pmf(self, k):
    """Probability mass function P(N = k)."""
    return eval_dist(lambda v: geom.pmf(v + 1, self.p), k)
cdf 
cdf(k)

Cumulative distribution function P(N <= k).

Source code in lossmodels/frequency/geometric.py 
46
47
48
def cdf(self, k):
    """Cumulative distribution function P(N <= k)."""
    return eval_dist(lambda v: geom.cdf(v + 1, self.p), k)

NegativeBinomial

Bases: FrequencyModel

Negative Binomial frequency model.

Parameterization

N = number of failures before the r-th success Support: {0, 1, 2, ...}

Parameters

r : float Number of successes, with r > 0. p : float Probability of success, with 0 < p <= 1.

Notes

This matches SciPy's negative binomial parameterization: scipy.stats.nbinom(r, p)

Under this convention

E[N] = r(1 - p) / p Var(N) = r(1 - p) / p^2

Source code in lossmodels/frequency/negbinomial.py 
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
class NegativeBinomial(FrequencyModel):
    """
    Negative Binomial frequency model.

    Parameterization
    ----------------
    N = number of failures before the r-th success
    Support: {0, 1, 2, ...}

    Parameters
    ----------
    r : float
        Number of successes, with r > 0.
    p : float
        Probability of success, with 0 < p <= 1.

    Notes
    -----
    This matches SciPy's negative binomial parameterization:
    scipy.stats.nbinom(r, p)

    Under this convention:
        E[N] = r(1 - p) / p
        Var(N) = r(1 - p) / p^2
    """

    def __init__(self, r: float, p: float):
        if r <= 0:
            raise ValueError("r must be positive.")
        if not (0 < p <= 1):
            raise ValueError("p must be in (0, 1].")

        self.r = r
        self.p = p

    def sample(self, size: int = 1) -> np.ndarray:
        """
        Generate random samples of claim counts.
        """
        if size <= 0:
            raise ValueError("size must be positive.")

        return np.random.negative_binomial(self.r, self.p, size=size)

    def mean(self) -> float:
        """
        E[N] = r(1 - p) / p
        """
        return self.r * (1 - self.p) / self.p

    def variance(self) -> float:
        """
        Var(N) = r(1 - p) / p^2
        """
        return self.r * (1 - self.p) / (self.p ** 2)

    def pmf(self, k):
        return eval_dist(lambda v: nbinom.pmf(v, self.r, self.p), k)

    def cdf(self, k):
        return eval_dist(lambda v: nbinom.cdf(v, self.r, self.p), k)

    def __repr__(self) -> str:
        return f"NegativeBinomial(r={self.r}, p={self.p})"
sample 
sample(size: int = 1) -> np.ndarray

Generate random samples of claim counts.

Source code in lossmodels/frequency/negbinomial.py 
43
44
45
46
47
48
49
50
def sample(self, size: int = 1) -> np.ndarray:
    """
    Generate random samples of claim counts.
    """
    if size <= 0:
        raise ValueError("size must be positive.")

    return np.random.negative_binomial(self.r, self.p, size=size)
mean 
mean() -> float

E[N] = r(1 - p) / p

Source code in lossmodels/frequency/negbinomial.py 
52
53
54
55
56
def mean(self) -> float:
    """
    E[N] = r(1 - p) / p
    """
    return self.r * (1 - self.p) / self.p
variance 
variance() -> float

Var(N) = r(1 - p) / p^2

Source code in lossmodels/frequency/negbinomial.py 
58
59
60
61
62
def variance(self) -> float:
    """
    Var(N) = r(1 - p) / p^2
    """
    return self.r * (1 - self.p) / (self.p ** 2)

Poisson

Bases: FrequencyModel

Poisson frequency model.

Parameterization

N ~ Poisson(lam)

Support: {0, 1, 2, ...}

Parameters

lam : float Expected claim count, with lam >= 0.

Source code in lossmodels/frequency/poisson.py 
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
class Poisson(FrequencyModel):
    """
    Poisson frequency model.

    Parameterization
    ----------------
    N ~ Poisson(lam)

    Support: {0, 1, 2, ...}

    Parameters
    ----------
    lam : float
        Expected claim count, with lam >= 0.
    """

    def __init__(self, lam: float):
        if lam < 0:
            raise ValueError("lam must be nonnegative.")
        self.lam = lam

    def sample(self, size: int = 1) -> np.ndarray:
        """Generate random samples of claim counts."""
        if size <= 0:
            raise ValueError("size must be positive.")
        return np.random.poisson(self.lam, size=size)

    def mean(self) -> float:
        """E[N] = lam."""
        return self.lam

    def variance(self) -> float:
        """Var(N) = lam."""
        return self.lam

    def pmf(self, k):
        """Probability mass function P(N = k)."""
        return eval_dist(lambda v: poisson.pmf(v, self.lam), k)

    def cdf(self, k):
        """Cumulative distribution function P(N <= k)."""
        return eval_dist(lambda v: poisson.cdf(v, self.lam), k)

    def __repr__(self) -> str:
        return f"Poisson(lam={self.lam})"
sample 
sample(size: int = 1) -> np.ndarray

Generate random samples of claim counts.

Source code in lossmodels/frequency/poisson.py 
29
30
31
32
33
def sample(self, size: int = 1) -> np.ndarray:
    """Generate random samples of claim counts."""
    if size <= 0:
        raise ValueError("size must be positive.")
    return np.random.poisson(self.lam, size=size)
mean 
mean() -> float

E[N] = lam.

Source code in lossmodels/frequency/poisson.py 
35
36
37
def mean(self) -> float:
    """E[N] = lam."""
    return self.lam
variance 
variance() -> float

Var(N) = lam.

Source code in lossmodels/frequency/poisson.py 
39
40
41
def variance(self) -> float:
    """Var(N) = lam."""
    return self.lam
pmf 
pmf(k)

Probability mass function P(N = k).

Source code in lossmodels/frequency/poisson.py 
43
44
45
def pmf(self, k):
    """Probability mass function P(N = k)."""
    return eval_dist(lambda v: poisson.pmf(v, self.lam), k)
cdf 
cdf(k)

Cumulative distribution function P(N <= k).

Source code in lossmodels/frequency/poisson.py 
47
48
49
def cdf(self, k):
    """Cumulative distribution function P(N <= k)."""
    return eval_dist(lambda v: poisson.cdf(v, self.lam), k)

ZeroTruncated

Bases: FrequencyModel

Zero-truncated version of an (a, b, 0) frequency model.

Parameters

base : FrequencyModel Any frequency model exposing pmf, cdf, mean, variance, and sample (e.g. Poisson(2.0)). Its mass at zero is removed and the remaining probabilities are rescaled by 1 / (1 - p_0).

Source code in lossmodels/frequency/truncated.py 
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
class ZeroTruncated(FrequencyModel):
    """Zero-truncated version of an (a, b, 0) frequency model.

    Parameters
    ----------
    base : FrequencyModel
        Any frequency model exposing ``pmf``, ``cdf``, ``mean``, ``variance``,
        and ``sample`` (e.g. ``Poisson(2.0)``). Its mass at zero is removed and
        the remaining probabilities are rescaled by ``1 / (1 - p_0)``.
    """

    def __init__(self, base: FrequencyModel):
        for attr in ("pmf", "cdf", "mean", "variance", "sample"):
            if not hasattr(base, attr):
                raise TypeError(f"base must implement {attr}().")
        self.base = base
        self._p0 = float(base.pmf(0))
        if self._p0 >= 1.0:
            raise ValueError("base places all mass at zero; cannot truncate.")

    def pmf(self, k):
        def f(v):
            v = np.asarray(v, dtype=float)
            return np.where(v >= 1, self.base.pmf(v) / (1.0 - self._p0), 0.0)

        return eval_dist(f, k)

    def cdf(self, k):
        def f(v):
            v = np.asarray(v, dtype=float)
            return np.where(
                v >= 1, (self.base.cdf(v) - self._p0) / (1.0 - self._p0), 0.0
            )

        return eval_dist(f, k)

    def mean(self) -> float:
        return self.base.mean() / (1.0 - self._p0)

    def variance(self) -> float:
        m = self.mean()
        e2 = (self.base.variance() + self.base.mean() ** 2) / (1.0 - self._p0)
        return e2 - m ** 2

    def sample(self, size: int = 1) -> np.ndarray:
        if size <= 0:
            raise ValueError("size must be positive.")
        pieces = []
        n = 0
        while n < size:
            draw = np.asarray(self.base.sample(max(2 * (size - n), 16)))
            nz = draw[draw > 0].astype(int)
            pieces.append(nz)
            n += nz.size
        return np.concatenate(pieces)[:size]

    def __repr__(self) -> str:
        return f"ZeroTruncated({self.base!r})"

ZeroModified

Bases: FrequencyModel

Zero-modified version of an (a, b, 0) frequency model.

Parameters

base : FrequencyModel Any frequency model exposing pmf, cdf, mean, variance, and sample. p0_modified : float The (arbitrary) probability mass placed at zero, in [0, 1). Setting it to 0 recovers the zero-truncated distribution.

Source code in lossmodels/frequency/modified.py 
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
class ZeroModified(FrequencyModel):
    """Zero-modified version of an (a, b, 0) frequency model.

    Parameters
    ----------
    base : FrequencyModel
        Any frequency model exposing ``pmf``, ``cdf``, ``mean``, ``variance``,
        and ``sample``.
    p0_modified : float
        The (arbitrary) probability mass placed at zero, in [0, 1). Setting it to
        0 recovers the zero-truncated distribution.
    """

    def __init__(self, base: FrequencyModel, p0_modified: float):
        for attr in ("pmf", "cdf", "mean", "variance", "sample"):
            if not hasattr(base, attr):
                raise TypeError(f"base must implement {attr}().")
        if not 0.0 <= p0_modified < 1.0:
            raise ValueError("p0_modified must be in [0, 1).")
        self.base = base
        self.p0_modified = float(p0_modified)
        self._p0 = float(base.pmf(0))
        if self._p0 >= 1.0:
            raise ValueError("base places all mass at zero; cannot modify.")

    def pmf(self, k):
        scale = (1.0 - self.p0_modified) / (1.0 - self._p0)

        def f(v):
            v = np.asarray(v, dtype=float)
            return np.where(
                v == 0,
                self.p0_modified,
                np.where(v >= 1, scale * self.base.pmf(v), 0.0),
            )

        return eval_dist(f, k)

    def cdf(self, k):
        scale = (1.0 - self.p0_modified) / (1.0 - self._p0)

        def f(v):
            v = np.asarray(v, dtype=float)
            upper = self.p0_modified + scale * (self.base.cdf(v) - self._p0)
            return np.where(v < 0, 0.0, np.where(v < 1, self.p0_modified, upper))

        return eval_dist(f, k)

    def _zt_moments(self):
        ztm = self.base.mean() / (1.0 - self._p0)
        ztv = (self.base.variance() + self.base.mean() ** 2) / (1.0 - self._p0) - ztm ** 2
        return ztm, ztv

    def mean(self) -> float:
        ztm, _ = self._zt_moments()
        return (1.0 - self.p0_modified) * ztm

    def variance(self) -> float:
        ztm, ztv = self._zt_moments()
        return (
            (1.0 - self.p0_modified) * ztv
            + self.p0_modified * (1.0 - self.p0_modified) * ztm ** 2
        )

    def sample(self, size: int = 1) -> np.ndarray:
        if size <= 0:
            raise ValueError("size must be positive.")
        keep = np.random.random(size) >= self.p0_modified
        n_keep = int(keep.sum())
        out = np.zeros(size, dtype=int)
        if n_keep:
            out[keep] = ZeroTruncated(self.base).sample(n_keep)
        return out

    def __repr__(self) -> str:
        return f"ZeroModified({self.base!r}, p0_modified={self.p0_modified})"

Logarithmic

Bases: FrequencyModel

Logarithmic frequency distribution (Loss Models B.3.1.3).

N ~ Logarithmic(beta), support {1, 2, ...}. p_k = (beta / (1+beta))^k / (k ln(1+beta)), k = 1, 2, ... E[N] = beta / ln(1+beta) Var[N] = beta [1 + beta - beta/ln(1+beta)] / ln(1+beta)

It is the r -> 0 limit of the zero-truncated negative binomial.

Source code in lossmodels/frequency/truncated.py 
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
class Logarithmic(FrequencyModel):
    """Logarithmic frequency distribution (Loss Models B.3.1.3).

    N ~ Logarithmic(beta), support {1, 2, ...}.
        p_k = (beta / (1+beta))^k / (k ln(1+beta)),  k = 1, 2, ...
        E[N] = beta / ln(1+beta)
        Var[N] = beta [1 + beta - beta/ln(1+beta)] / ln(1+beta)

    It is the r -> 0 limit of the zero-truncated negative binomial.
    """

    def __init__(self, beta: float):
        if beta <= 0:
            raise ValueError("beta must be positive.")
        self.beta = beta
        self._c = log(1.0 + beta)
        self._r = beta / (1.0 + beta)

    def pmf(self, k):
        def f(v):
            v = np.asarray(v, dtype=float)
            is_pos_int = (v >= 1) & (np.mod(v, 1) == 0)
            safe = np.where(is_pos_int, v, 1.0)
            return np.where(is_pos_int, self._r ** safe / (safe * self._c), 0.0)

        return eval_dist(f, k)

    def cdf(self, k):
        def f(v):
            arr = np.asarray(v, dtype=float)
            flat = np.atleast_1d(arr)
            idx = np.floor(flat).astype(int)
            kmax = int(idx.max()) if idx.size and idx.max() >= 1 else 0
            if kmax >= 1:
                js = np.arange(1, kmax + 1, dtype=float)
                cum = np.concatenate([[0.0], np.cumsum(self._r ** js / (js * self._c))])
            else:
                cum = np.array([0.0])
            out = np.where(idx < 1, 0.0, cum[np.clip(idx, 0, kmax)])
            return out.reshape(arr.shape)

        return eval_dist(f, k)

    def mean(self) -> float:
        return self.beta / self._c

    def variance(self) -> float:
        return self.beta * (1.0 + self.beta - self.beta / self._c) / self._c

    def sample(self, size: int = 1) -> np.ndarray:
        if size <= 0:
            raise ValueError("size must be positive.")
        return np.random.logseries(self._r, size=size)

    def __repr__(self) -> str:
        return f"Logarithmic(beta={self.beta})"

CollectiveRiskModel

Bases: AggregateModel

Collective risk model for aggregate loss:

S = X1 + X2 + ... + XN

where: - N is the claim count random variable (frequency) - Xi are iid claim severities

Assumes: - severities are iid - N is independent of severities

Source code in lossmodels/aggregate/collective.py 
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
class CollectiveRiskModel(AggregateModel):
    """
    Collective risk model for aggregate loss:

        S = X1 + X2 + ... + XN

    where:
    - N is the claim count random variable (frequency)
    - Xi are iid claim severities

    Assumes:
    - severities are iid
    - N is independent of severities
    """

    def __init__(self, frequency: FrequencyModel, severity: SeverityModel):
        self.frequency = frequency
        self.severity = severity

    def sample(self, size: int = 1) -> np.ndarray:
        """
        Generate random samples of aggregate loss.
        """
        if size <= 0:
            raise ValueError("size must be positive")

        counts = self.frequency.sample(size=size)
        aggregate_losses = np.zeros(size, dtype=float)

        for i, n_claims in enumerate(counts):
            if n_claims > 0:
                aggregate_losses[i] = np.sum(self.severity.sample(size=int(n_claims)))

        return aggregate_losses

    def mean(self) -> float:
        """
        E[S] = E[N] * E[X]
        """
        return self.frequency.mean() * self.severity.mean()

    def variance(self) -> float:
        """
        Var(S) = E[N] Var(X) + Var(N) (E[X])^2
        """
        en = self.frequency.mean()
        vn = self.frequency.variance()
        ex = self.severity.mean()
        vx = self.severity.variance()

        return en * vx + vn * (ex ** 2)

    def frequency_mean(self) -> float:
        return self.frequency.mean()

    def severity_mean(self) -> float:
        return self.severity.mean()

    def summary(self) -> dict:
        """
        Return a small summary of the model.
        """
        return {
            "frequency_model": repr(self.frequency),
            "severity_model": repr(self.severity),
            "frequency_mean": self.frequency.mean(),
            "severity_mean": self.severity.mean(),
            "aggregate_mean": self.mean(),
            "aggregate_variance": self.variance(),
            "aggregate_std": self.std(),
        }

    def __repr__(self):
        return (
            f"CollectiveRiskModel("
            f"frequency={repr(self.frequency)}, "
            f"severity={repr(self.severity)})"
        )
sample 
sample(size: int = 1) -> np.ndarray

Generate random samples of aggregate loss.

Source code in lossmodels/aggregate/collective.py 
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
def sample(self, size: int = 1) -> np.ndarray:
    """
    Generate random samples of aggregate loss.
    """
    if size <= 0:
        raise ValueError("size must be positive")

    counts = self.frequency.sample(size=size)
    aggregate_losses = np.zeros(size, dtype=float)

    for i, n_claims in enumerate(counts):
        if n_claims > 0:
            aggregate_losses[i] = np.sum(self.severity.sample(size=int(n_claims)))

    return aggregate_losses
mean 
mean() -> float

E[S] = E[N] * E[X]

Source code in lossmodels/aggregate/collective.py 
43
44
45
46
47
def mean(self) -> float:
    """
    E[S] = E[N] * E[X]
    """
    return self.frequency.mean() * self.severity.mean()
variance 
variance() -> float

Var(S) = E[N] Var(X) + Var(N) (E[X])^2

Source code in lossmodels/aggregate/collective.py 
49
50
51
52
53
54
55
56
57
58
def variance(self) -> float:
    """
    Var(S) = E[N] Var(X) + Var(N) (E[X])^2
    """
    en = self.frequency.mean()
    vn = self.frequency.variance()
    ex = self.severity.mean()
    vx = self.severity.variance()

    return en * vx + vn * (ex ** 2)
summary 
summary() -> dict

Return a small summary of the model.

Source code in lossmodels/aggregate/collective.py 
66
67
68
69
70
71
72
73
74
75
76
77
78
def summary(self) -> dict:
    """
    Return a small summary of the model.
    """
    return {
        "frequency_model": repr(self.frequency),
        "severity_model": repr(self.severity),
        "frequency_mean": self.frequency.mean(),
        "severity_mean": self.severity.mean(),
        "aggregate_mean": self.mean(),
        "aggregate_variance": self.variance(),
        "aggregate_std": self.std(),
    }

EmpiricalSeverity

Bases: SeverityModel

Empirical severity model based on observed loss data.

Parameters

data : array-like Observed severity values. Must be nonempty and nonnegative.

Source code in lossmodels/empirical/distribution.py 
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
class EmpiricalSeverity(SeverityModel):
    """
    Empirical severity model based on observed loss data.

    Parameters
    ----------
    data : array-like
        Observed severity values. Must be nonempty and nonnegative.
    """

    def __init__(self, data):
        data = np.asarray(data, dtype=float)

        if data.size == 0:
            raise ValueError("data must not be empty.")
        if np.any(data < 0):
            raise ValueError("severity data must be nonnegative.")

        self.data = data

    def sample(self, size: int = 1) -> np.ndarray:
        """
        Generate bootstrap samples from the empirical severity distribution.
        """
        if size <= 0:
            raise ValueError("size must be positive.")

        return np.random.choice(self.data, size=size, replace=True)

    def mean(self) -> float:
        return float(np.mean(self.data))

    def variance(self) -> float:
        return float(np.var(self.data, ddof=0))

    def pdf(self, x: float) -> float:
        """
        Empirical severity is discrete, so this returns the empirical point mass at x.
        For continuous-looking data, this will often be 0 except at exact observed values.
        """
        if x < 0:
            return 0.0

        return float(np.mean(self.data == x))

    def cdf(self, x: float) -> float:
        if x < 0:
            return 0.0

        return float(np.mean(self.data <= x))

    def quantile(self, p):
        """Empirical quantile (inverse CDF) from the observed data."""
        q = np.quantile(self.data, p)
        return float(q) if np.ndim(p) == 0 else np.asarray(q, dtype=float)

    def excess_loss(self, d: float) -> float:
        """
        E[(X - d)+] computed empirically.
        """
        if d < 0:
            raise ValueError("d must be nonnegative.")

        return float(np.mean(np.maximum(self.data - d, 0.0)))

    def limited_expected_value(self, d: float) -> float:
        """
        E[min(X, d)] computed empirically.
        """
        if d < 0:
            raise ValueError("d must be nonnegative.")

        return float(np.mean(np.minimum(self.data, d)))

    def __repr__(self) -> str:
        return f"EmpiricalSeverity(n={len(self.data)})"
sample 
sample(size: int = 1) -> np.ndarray

Generate bootstrap samples from the empirical severity distribution.

Source code in lossmodels/empirical/distribution.py 
27
28
29
30
31
32
33
34
def sample(self, size: int = 1) -> np.ndarray:
    """
    Generate bootstrap samples from the empirical severity distribution.
    """
    if size <= 0:
        raise ValueError("size must be positive.")

    return np.random.choice(self.data, size=size, replace=True)
pdf 
pdf(x: float) -> float

Empirical severity is discrete, so this returns the empirical point mass at x. For continuous-looking data, this will often be 0 except at exact observed values.

Source code in lossmodels/empirical/distribution.py 
42
43
44
45
46
47
48
49
50
def pdf(self, x: float) -> float:
    """
    Empirical severity is discrete, so this returns the empirical point mass at x.
    For continuous-looking data, this will often be 0 except at exact observed values.
    """
    if x < 0:
        return 0.0

    return float(np.mean(self.data == x))
quantile 
quantile(p)

Empirical quantile (inverse CDF) from the observed data.

Source code in lossmodels/empirical/distribution.py 
58
59
60
61
def quantile(self, p):
    """Empirical quantile (inverse CDF) from the observed data."""
    q = np.quantile(self.data, p)
    return float(q) if np.ndim(p) == 0 else np.asarray(q, dtype=float)
excess_loss 
excess_loss(d: float) -> float

E[(X - d)+] computed empirically.

Source code in lossmodels/empirical/distribution.py 
63
64
65
66
67
68
69
70
def excess_loss(self, d: float) -> float:
    """
    E[(X - d)+] computed empirically.
    """
    if d < 0:
        raise ValueError("d must be nonnegative.")

    return float(np.mean(np.maximum(self.data - d, 0.0)))
limited_expected_value 
limited_expected_value(d: float) -> float

E[min(X, d)] computed empirically.

Source code in lossmodels/empirical/distribution.py 
72
73
74
75
76
77
78
79
def limited_expected_value(self, d: float) -> float:
    """
    E[min(X, d)] computed empirically.
    """
    if d < 0:
        raise ValueError("d must be nonnegative.")

    return float(np.mean(np.minimum(self.data, d)))

EmpiricalFrequency

Bases: FrequencyModel

Empirical frequency model based on observed claim count data.

Parameters

data : array-like Observed claim counts. Must be nonempty, nonnegative, and integer-valued.

Source code in lossmodels/empirical/distribution.py 
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
class EmpiricalFrequency(FrequencyModel):
    """
    Empirical frequency model based on observed claim count data.

    Parameters
    ----------
    data : array-like
        Observed claim counts. Must be nonempty, nonnegative, and integer-valued.
    """

    def __init__(self, data):
        data = np.asarray(data)

        if data.size == 0:
            raise ValueError("data must not be empty.")
        if np.any(data < 0):
            raise ValueError("frequency data must be nonnegative.")
        if not np.all(np.equal(data, np.floor(data))):
            raise ValueError("frequency data must be integer-valued.")

        self.data = data.astype(int)

    def sample(self, size: int = 1) -> np.ndarray:
        """
        Generate bootstrap samples from the empirical frequency distribution.
        """
        if size <= 0:
            raise ValueError("size must be positive.")

        return np.random.choice(self.data, size=size, replace=True)

    def mean(self) -> float:
        return float(np.mean(self.data))

    def variance(self) -> float:
        return float(np.var(self.data, ddof=0))

    def pmf(self, k: int) -> float:
        if k < 0:
            return 0.0

        return float(np.mean(self.data == k))

    def cdf(self, k: int) -> float:
        if k < 0:
            return 0.0

        return float(np.mean(self.data <= k))

    def __repr__(self) -> str:
        return f"EmpiricalFrequency(n={len(self.data)})"
sample 
sample(size: int = 1) -> np.ndarray

Generate bootstrap samples from the empirical frequency distribution.

Source code in lossmodels/empirical/distribution.py 
107
108
109
110
111
112
113
114
def sample(self, size: int = 1) -> np.ndarray:
    """
    Generate bootstrap samples from the empirical frequency distribution.
    """
    if size <= 0:
        raise ValueError("size must be positive.")

    return np.random.choice(self.data, size=size, replace=True)

fit_exponential 

fit_exponential(data) -> Exponential

Fit an Exponential severity model by maximum likelihood.

For X_i ~ Exponential(rate), the MLE is: rate_hat = 1 / mean(data)

Source code in lossmodels/estimation/mle.py 
36
37
38
39
40
41
42
43
44
45
46
47
48
def fit_exponential(data) -> Exponential:
    """
    Fit an Exponential severity model by maximum likelihood.

    For X_i ~ Exponential(rate), the MLE is:
        rate_hat = 1 / mean(data)
    """
    data = _validate_positive_data(data)
    mean_x = float(np.mean(data))
    if mean_x <= 0:
        raise ValueError("Mean of data must be positive.")
    rate_hat = 1.0 / mean_x
    return Exponential(rate=float(rate_hat))

fit_gamma 

fit_gamma(data) -> Gamma

Fit a Gamma severity model by maximum likelihood using SciPy.

Returns

Gamma Fitted Gamma(alpha, theta) model.

Notes

This constrains loc = 0 so the support is x > 0, consistent with the severity model implementation.

Source code in lossmodels/estimation/mle.py 
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
def fit_gamma(data) -> Gamma:
    """
    Fit a Gamma severity model by maximum likelihood using SciPy.

    Returns
    -------
    Gamma
        Fitted Gamma(alpha, theta) model.

    Notes
    -----
    This constrains loc = 0 so the support is x > 0, consistent with the
    severity model implementation.
    """
    data = _validate_positive_data(data)
    alpha_hat, loc_hat, theta_hat = gamma_dist.fit(data, floc=0)
    if loc_hat != 0:
        raise RuntimeError("Gamma fit returned nonzero location despite floc=0.")
    return Gamma(alpha=float(alpha_hat), theta=float(theta_hat))

fit_lognormal 

fit_lognormal(data) -> Lognormal

Fit a Lognormal severity model by maximum likelihood.

If log(X) ~ Normal(mu, sigma^2), the MLEs are: mu_hat = mean(log(data)) sigma_hat = sqrt(mean((log(data) - mu_hat)^2))

Notes

This uses the MLE version of the variance (ddof=0).

Source code in lossmodels/estimation/mle.py 
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
def fit_lognormal(data) -> Lognormal:
    """
    Fit a Lognormal severity model by maximum likelihood.

    If log(X) ~ Normal(mu, sigma^2), the MLEs are:
        mu_hat = mean(log(data))
        sigma_hat = sqrt(mean((log(data) - mu_hat)^2))

    Notes
    -----
    This uses the MLE version of the variance (ddof=0).
    """
    data = _validate_positive_data(data)
    log_data = np.log(data)
    mu_hat = float(np.mean(log_data))
    sigma_hat = float(np.sqrt(np.mean((log_data - mu_hat) ** 2)))
    return Lognormal(mu=mu_hat, sigma=sigma_hat)

fit_pareto 

fit_pareto(data) -> Pareto

Fit a Pareto Type I severity model by maximum likelihood.

For X_i ~ Pareto(alpha, theta) with support x >= theta, the MLEs are:

theta_hat = min(data)
alpha_hat = n / sum(log(data / theta_hat))
Returns

Pareto Fitted Pareto(alpha, theta) model.

Source code in lossmodels/estimation/mle.py 
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
def fit_pareto(data) -> Pareto:
    """
    Fit a Pareto Type I severity model by maximum likelihood.

    For X_i ~ Pareto(alpha, theta) with support x >= theta, the MLEs are:

        theta_hat = min(data)
        alpha_hat = n / sum(log(data / theta_hat))

    Returns
    -------
    Pareto
        Fitted Pareto(alpha, theta) model.
    """
    data = _validate_positive_data(data)

    theta_hat = float(np.min(data))
    if theta_hat <= 0:
        raise ValueError("Estimated theta must be positive.")

    log_ratios = np.log(data / theta_hat)
    denom = float(np.sum(log_ratios))
    if denom <= 0:
        raise ValueError(
            "Pareto MLE requires data with at least one observation above the minimum."
        )

    alpha_hat = float(len(data) / denom)
    return Pareto(alpha=alpha_hat, theta=theta_hat)

fit_poisson 

fit_poisson(data) -> Poisson

Fit a Poisson frequency model by maximum likelihood.

For N_i ~ Poisson(lam), the MLE is: lam_hat = mean(data)

Notes

An all-zero dataset is valid and yields lam_hat = 0.

Source code in lossmodels/estimation/mle.py 
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
def fit_poisson(data) -> Poisson:
    """
    Fit a Poisson frequency model by maximum likelihood.

    For N_i ~ Poisson(lam), the MLE is:
        lam_hat = mean(data)

    Notes
    -----
    An all-zero dataset is valid and yields lam_hat = 0.
    """
    data = _validate_count_data(data)
    lam_hat = float(np.mean(data))
    if lam_hat < 0:
        raise ValueError("Estimated lambda must be nonnegative.")
    return Poisson(lam=lam_hat)

fit_weibull 

fit_weibull(data) -> Weibull

Fit a Weibull severity model by maximum likelihood using SciPy.

Returns

Weibull Fitted Weibull(k, lam) model.

Notes

This constrains loc = 0 so the support is x > 0, consistent with the severity model implementation.

Source code in lossmodels/estimation/mle.py 
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
def fit_weibull(data) -> Weibull:
    """
    Fit a Weibull severity model by maximum likelihood using SciPy.

    Returns
    -------
    Weibull
        Fitted Weibull(k, lam) model.

    Notes
    -----
    This constrains loc = 0 so the support is x > 0, consistent with the
    severity model implementation.
    """
    data = _validate_positive_data(data)
    k_hat, loc_hat, lam_hat = weibull_min.fit(data, floc=0)
    if loc_hat != 0:
        raise RuntimeError("Weibull fit returned nonzero location despite floc=0.")
    return Weibull(k=float(k_hat), lam=float(lam_hat))

fit_negbinomial 

fit_negbinomial(data) -> NegativeBinomial

Fit a Negative Binomial frequency model by numerical maximum likelihood.

Parameterization

N = number of failures before the r-th success Support: {0, 1, 2, ...} Mean = r(1-p)/p Variance = r(1-p)/p^2

Source code in lossmodels/estimation/mle.py 
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
def fit_negbinomial(data) -> NegativeBinomial:
    """
    Fit a Negative Binomial frequency model by numerical maximum likelihood.

    Parameterization
    ----------------
    N = number of failures before the r-th success
    Support: {0, 1, 2, ...}
    Mean = r(1-p)/p
    Variance = r(1-p)/p^2
    """
    data = _validate_count_data(data)

    mean_x = float(np.mean(data))
    var_x = float(np.var(data, ddof=0))

    if var_x > mean_x and mean_x > 0:
        p0 = mean_x / var_x
        r0 = mean_x**2 / (var_x - mean_x)
        initial = np.array([r0, p0], dtype=float)
    else:
        initial = np.array([1.0, 0.5], dtype=float)

    bounds = [
        (1e-8, None),           # r > 0
        (1e-8, 1.0 - 1e-8),     # 0 < p < 1
    ]

    def neg_log_likelihood(params):
        r, p = params
        try:
            model = NegativeBinomial(r=r, p=p)
            pmf_vals = np.array([model.pmf(int(x)) for x in data], dtype=float)
            if np.any(~np.isfinite(pmf_vals)) or np.any(pmf_vals <= 0):
                return np.inf
            return float(-np.sum(np.log(pmf_vals)))
        except Exception:
            return np.inf

    result = minimize(
        neg_log_likelihood,
        x0=initial,
        bounds=bounds,
        method="L-BFGS-B",
    )

    if not result.success:
        raise RuntimeError(
            f"Negative Binomial MLE optimization failed: {result.message}"
        )

    r_hat, p_hat = result.x
    return NegativeBinomial(r=float(r_hat), p=float(p_hat))

fit_mle 

fit_mle(model_class, data, initial_params, bounds=None)

Generic numerical maximum likelihood estimation for models with a pdf method.

Parameters

model_class : class A model class that can be instantiated as model_class(*params) and provides a pdf(x) method. data : array-like Observed data. initial_params : array-like Initial parameter guess for the optimizer. bounds : list of tuple, optional Bounds passed to scipy.optimize.minimize.

Returns

object Fitted model instance of type model_class.

Source code in lossmodels/estimation/mle.py 
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
def fit_mle(model_class, data, initial_params, bounds=None):
    """
    Generic numerical maximum likelihood estimation for models with a pdf method.

    Parameters
    ----------
    model_class : class
        A model class that can be instantiated as model_class(*params)
        and provides a pdf(x) method.
    data : array-like
        Observed data.
    initial_params : array-like
        Initial parameter guess for the optimizer.
    bounds : list of tuple, optional
        Bounds passed to scipy.optimize.minimize.

    Returns
    -------
    object
        Fitted model instance of type model_class.
    """
    data = _validate_positive_data(data)
    initial_params = np.asarray(initial_params, dtype=float)

    if initial_params.size == 0:
        raise ValueError("initial_params must not be empty.")

    def neg_log_likelihood(params):
        try:
            model = model_class(*params)
            pdf_vals = np.array([model.pdf(x) for x in data], dtype=float)
            if np.any(~np.isfinite(pdf_vals)) or np.any(pdf_vals <= 0):
                return np.inf
            return float(-np.sum(np.log(pdf_vals)))
        except Exception:
            return np.inf

    result = minimize(
        neg_log_likelihood,
        x0=initial_params,
        bounds=bounds,
        method="L-BFGS-B" if bounds is not None else "BFGS",
    )

    if not result.success:
        raise RuntimeError(f"MLE optimization failed: {result.message}")

    return model_class(*result.x)

fit_best_severity 

fit_best_severity(
    data, candidates=None, method="mle", criterion="aic"
)

Fit a set of severity models and return the best one by AIC or BIC.

Parameters

data : array-like Severity observations. candidates : list of str, optional Candidate model names. Defaults to all supported severity fitters. method : {"mle", "moments"} Fitting method. criterion : {"aic", "bic"} Selection criterion.

Returns

dict Dictionary with keys: - "best_name" - "best_model" - "criterion" - "method" - "results"

Notes

Models that fail to fit are skipped.

Source code in lossmodels/estimation/model_selection.py 
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
def fit_best_severity(data, candidates=None, method="mle", criterion="aic"):
    """
    Fit a set of severity models and return the best one by AIC or BIC.

    Parameters
    ----------
    data : array-like
        Severity observations.
    candidates : list of str, optional
        Candidate model names. Defaults to all supported severity fitters.
    method : {"mle", "moments"}
        Fitting method.
    criterion : {"aic", "bic"}
        Selection criterion.

    Returns
    -------
    dict
        Dictionary with keys:
        - "best_name"
        - "best_model"
        - "criterion"
        - "method"
        - "results"

    Notes
    -----
    Models that fail to fit are skipped.
    """
    data = np.asarray(data, dtype=float)
    if data.size == 0:
        raise ValueError("data must not be empty.")

    if method not in {"mle", "moments"}:
        raise ValueError("method must be 'mle' or 'moments'.")
    if criterion not in {"aic", "bic"}:
        raise ValueError("criterion must be 'aic' or 'bic'.")

    fitters = SEVERITY_MLE_FITTERS if method == "mle" else SEVERITY_MOMENT_FITTERS

    if candidates is None:
        candidates = list(fitters.keys())

    results = []

    for name in candidates:
        if name not in fitters:
            raise ValueError(f"Unsupported candidate: {name}")

        fitter, k = fitters[name]

        try:
            model = fitter(data)
            score = aic(model, data, k=k) if criterion == "aic" else bic(model, data, k=k)

            results.append({
                "name": name,
                "model": model,
                "k": k,
                "score": float(score),
            })
        except Exception as exc:
            results.append({
                "name": name,
                "model": None,
                "k": k,
                "score": float("inf"),
                "error": str(exc),
            })

    valid = [r for r in results if np.isfinite(r["score"])]
    if not valid:
        error_details = "; ".join(
            f'{r["name"]}: {r.get("error", "unknown error")}'
            for r in results
        )
        raise RuntimeError(
            f"No candidate models could be fit successfully. Errors: {error_details}"
        )

    best = min(valid, key=lambda r: r["score"])

    return {
        "best_name": best["name"],
        "best_model": best["model"],
        "criterion": criterion,
        "method": method,
        "results": results,
    }

fit_best_frequency 

fit_best_frequency(
    data, candidates=None, method="mle", criterion="aic"
)

Fit a set of frequency models and return the best one by AIC or BIC.

Source code in lossmodels/estimation/frequency_selection.py 
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
def fit_best_frequency(data, candidates=None, method="mle", criterion="aic"):
    """
    Fit a set of frequency models and return the best one by AIC or BIC.
    """
    data = np.asarray(data)

    if data.size == 0:
        raise ValueError("data must not be empty.")
    if np.any(data < 0):
        raise ValueError("frequency data must be nonnegative.")

    if method not in {"mle", "moments"}:
        raise ValueError("method must be 'mle' or 'moments'.")
    if criterion not in {"aic", "bic"}:
        raise ValueError("criterion must be 'aic' or 'bic'.")

    fitters = (
        FREQUENCY_MLE_FITTERS if method == "mle" else FREQUENCY_MOMENT_FITTERS
    )

    if candidates is None:
        candidates = list(fitters.keys())

    results = []

    for name in candidates:
        if name not in fitters:
            raise ValueError(f"Unsupported candidate: {name}")

        fitter, k = fitters[name]

        try:
            model = fitter(data)
            score = aic(model, data, k=k) if criterion == "aic" else bic(model, data, k=k)

            results.append({
                "name": name,
                "model": model,
                "k": k,
                "score": float(score),
            })

        except Exception as exc:
            results.append({
                "name": name,
                "model": None,
                "k": k,
                "score": float("inf"),
                "error": str(exc),
            })

    valid = [r for r in results if np.isfinite(r["score"])]
    if not valid:
        raise RuntimeError("No candidate frequency models could be fit successfully.")

    best = min(valid, key=lambda r: r["score"])

    return {
        "best_name": best["name"],
        "best_model": best["model"],
        "criterion": criterion,
        "method": method,
        "results": results,
    }

log_likelihood 

log_likelihood(model, data) -> float

Compute the log-likelihood of observed data under a fitted model.

Parameters

model : object Model instance with a pdf(x) method for continuous models or pmf(x) method for discrete models. data : array-like Observed data.

Returns

float Log-likelihood value.

Source code in lossmodels/estimation/diagnostics.py 
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
def log_likelihood(model, data) -> float:
    """
    Compute the log-likelihood of observed data under a fitted model.

    Parameters
    ----------
    model : object
        Model instance with a pdf(x) method for continuous models
        or pmf(x) method for discrete models.
    data : array-like
        Observed data.

    Returns
    -------
    float
        Log-likelihood value.
    """
    data = np.asarray(data)
    if data.size == 0:
        raise ValueError("data must not be empty.")

    if hasattr(model, "pdf"):
        evaluate = model.pdf
    elif hasattr(model, "pmf"):
        evaluate = model.pmf
    else:
        raise TypeError("model must implement either pdf(x) or pmf(x).")

    # Fast vectorized path; fall back to per-element for non-vectorized models.
    try:
        vals = np.asarray(evaluate(data), dtype=float)
        if vals.shape != data.shape:
            raise ValueError("non-vectorized result")
    except Exception:
        vals = np.array([evaluate(float(x)) for x in data], dtype=float)

    if np.any(~np.isfinite(vals)) or np.any(vals <= 0):
        return float(-np.inf)

    return float(np.sum(np.log(vals)))

aic 

aic(model, data, k: int) -> float

Compute Akaike Information Criterion.

Parameters

model : object Fitted model. data : array-like Observed data. k : int Number of estimated parameters.

Returns

float AIC value.

Source code in lossmodels/estimation/diagnostics.py 
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
def aic(model, data, k: int) -> float:
    """
    Compute Akaike Information Criterion.

    Parameters
    ----------
    model : object
        Fitted model.
    data : array-like
        Observed data.
    k : int
        Number of estimated parameters.

    Returns
    -------
    float
        AIC value.
    """
    if k <= 0:
        raise ValueError("k must be positive.")

    ll = log_likelihood(model, data)
    if not np.isfinite(ll):
        return float(np.inf)

    return float(2 * k - 2 * ll)

bic 

bic(model, data, k: int) -> float

Compute Bayesian Information Criterion.

Parameters

model : object Fitted model. data : array-like Observed data. k : int Number of estimated parameters.

Returns

float BIC value.

Source code in lossmodels/estimation/diagnostics.py 
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
def bic(model, data, k: int) -> float:
    """
    Compute Bayesian Information Criterion.

    Parameters
    ----------
    model : object
        Fitted model.
    data : array-like
        Observed data.
    k : int
        Number of estimated parameters.

    Returns
    -------
    float
        BIC value.
    """
    data = np.asarray(data)
    if data.size == 0:
        raise ValueError("data must not be empty.")
    if k <= 0:
        raise ValueError("k must be positive.")

    ll = log_likelihood(model, data)
    if not np.isfinite(ll):
        return float(np.inf)

    n = data.size
    return float(np.log(n) * k - 2 * ll)

ks_statistic 

ks_statistic(model, data) -> float

Kolmogorov-Smirnov distance sup_x |F_n(x) - F(x)| (smaller is better).

Most sensitive in the body of the distribution. See the module note on estimated-parameter p-values.

Source code in lossmodels/estimation/diagnostics.py 
138
139
140
141
142
143
144
145
146
147
148
149
def ks_statistic(model, data) -> float:
    """Kolmogorov-Smirnov distance sup_x |F_n(x) - F(x)| (smaller is better).

    Most sensitive in the body of the distribution. See the module note on
    estimated-parameter p-values.
    """
    x = _sorted_data(data)
    n = x.size
    f = _eval_cdf(model, x)
    d_plus = np.max(np.arange(1, n + 1) / n - f)
    d_minus = np.max(f - np.arange(0, n) / n)
    return float(max(d_plus, d_minus))

anderson_darling 

anderson_darling(model, data) -> float

Anderson-Darling statistic A^2 (smaller is better).

Weights the tails more than KS, so it is the better whole-distribution statistic for heavy-tailed loss data.

Source code in lossmodels/estimation/diagnostics.py 
152
153
154
155
156
157
158
159
160
161
162
163
def anderson_darling(model, data) -> float:
    """Anderson-Darling statistic A^2 (smaller is better).

    Weights the tails more than KS, so it is the better whole-distribution
    statistic for heavy-tailed loss data.
    """
    x = _sorted_data(data)
    n = x.size
    f = np.clip(_eval_cdf(model, x), 1e-12, 1.0 - 1e-12)
    i = np.arange(1, n + 1)
    s = np.sum((2 * i - 1) * (np.log(f) + np.log(1.0 - f[::-1])))
    return float(-n - s / n)

cramer_von_mises 

cramer_von_mises(model, data) -> float

Cramer-von Mises statistic W^2 (smaller is better).

Source code in lossmodels/estimation/diagnostics.py 
166
167
168
169
170
171
172
def cramer_von_mises(model, data) -> float:
    """Cramer-von Mises statistic W^2 (smaller is better)."""
    x = _sorted_data(data)
    n = x.size
    f = _eval_cdf(model, x)
    i = np.arange(1, n + 1)
    return float(1.0 / (12.0 * n) + np.sum((f - (2 * i - 1) / (2.0 * n)) ** 2))

tail_quantile_table 

tail_quantile_table(
    model, data, probs=(0.9, 0.95, 0.99, 0.995)
) -> list

Compare fitted vs empirical high quantiles -- the tail-fit check.

Returns a list of dicts with keys prob, empirical, fitted, abs_error, rel_error. A model can win on AIC and still miss here; for tail risk (VaR / TVaR / stop-loss) this is the diagnostic that matters. Requires the model to implement quantile(p).

Source code in lossmodels/estimation/diagnostics.py 
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
def tail_quantile_table(model, data, probs=(0.90, 0.95, 0.99, 0.995)) -> list:
    """Compare fitted vs empirical high quantiles -- the tail-fit check.

    Returns a list of dicts with keys ``prob``, ``empirical``, ``fitted``,
    ``abs_error``, ``rel_error``. A model can win on AIC and still miss here;
    for tail risk (VaR / TVaR / stop-loss) this is the diagnostic that matters.
    Requires the model to implement ``quantile(p)``.
    """
    if not hasattr(model, "quantile"):
        raise TypeError("model must implement quantile(p) for the tail table.")
    data = np.asarray(data, dtype=float)
    if data.size == 0:
        raise ValueError("data must not be empty.")
    rows = []
    for p in probs:
        emp = float(np.quantile(data, p))
        fit = float(model.quantile(p))
        rows.append({
            "prob": float(p),
            "empirical": emp,
            "fitted": fit,
            "abs_error": fit - emp,
            "rel_error": (fit - emp) / emp if emp != 0 else float("nan"),
        })
    return rows

goodness_of_fit 

goodness_of_fit(model, data, k: int) -> dict

One-call fit report combining relative and absolute measures.

Returns n, log_likelihood, aic, bic (relative -- compare across candidates) and ks, anderson_darling, cramer_von_mises (absolute distance-to-empirical -- smaller is better). See the module note on the estimated-parameter caveat for the distance statistics.

Source code in lossmodels/estimation/diagnostics.py 
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
def goodness_of_fit(model, data, k: int) -> dict:
    """One-call fit report combining relative and absolute measures.

    Returns ``n``, ``log_likelihood``, ``aic``, ``bic`` (relative -- compare
    across candidates) and ``ks``, ``anderson_darling``, ``cramer_von_mises``
    (absolute distance-to-empirical -- smaller is better). See the module note
    on the estimated-parameter caveat for the distance statistics.
    """
    return {
        "n": int(np.asarray(data).size),
        "log_likelihood": log_likelihood(model, data),
        "aic": aic(model, data, k=k),
        "bic": bic(model, data, k=k),
        "ks": ks_statistic(model, data),
        "anderson_darling": anderson_darling(model, data),
        "cramer_von_mises": cramer_von_mises(model, data),
    }