Density and intensity¶

fdatools.density: penalised log-density (fit_density) and log-intensity (fit_intensity) estimation with a basis expansion for the log function, fitted by damped Newton with exact integrals.

Penalised density and intensity estimation.

Two estimators replace R fda's density.fd (dropped from fda 6.3.0) and intensity.fd. Both write a positive function as the exponential of a basis expansion W(x) = Σ_k c_k φ_k(x) and maximise a penalised log-likelihood.

:func:fit_density estimates a probability density p(x) = C exp W(x) from a sample x_1 .. x_n by minimising

.. math::

F(c) = -\sum_{i=1}^{n} W(x_i) + n \log \int e^{W(u)}\,du
       + \lambda \int (LW)^2 ,

with C = 1 / ∫ e^W the normalising constant (Ramsay, Hooker & Graves 2009, § 5.4; Silverman 1986 for penalised log-density estimation).

:func:fit_intensity estimates the intensity μ(t) = exp W(t) of a non-homogeneous Poisson process from its event times by minimising the penalised negative log-likelihood of the process,

.. math::

F(c) = -\sum_{i=1}^{n} W(x_i) + \int e^{W(u)}\,du + \lambda \int (LW)^2 .

The penalty is λ cᵀRc with R the roughness matrix of the linear differential operator L (:meth:fdatools.Basis.penalty), the same scaling as R. Both criteria are convex in c. They are minimised by Newton's method with the exact Hessian and a backtracking (Armijo) line search; the integrals are computed by composite Gauss-Legendre quadrature, 12 nodes on each of the basis's refined quadrature panels, which is exact to rounding for the exponential of a polynomial piece of the size met in practice.

Convergence: the iteration stops once a full Newton step, which is then taken, moves no coefficient by more than tol * (1 + max|c|) (default tol = 1e-10); Newton's quadratic convergence leaves the result about tol² from the optimum. If rounding hides every further decrease of F first, the fit counts as converged when the last proposed step was within sqrt(tol) * (1 + max|c|).

When the basis reproduces the constant function and the penalty does not charge it (every derivative operator Dᵐ with m >= 1), the density criterion does not change under W -> W + s. The fit then returns the member of that family with ∫ W = 0 over the domain; the density itself, C exp W, is unique.

Examples:

>>> import numpy as np
>>> import fdatools as fdt
>>> from fdatools.density import fit_density
>>> rng = np.random.default_rng(0)
>>> x = rng.normal(size=300)
>>> basis = fdt.BSpline(domain=(-5.0, 5.0), n_basis=11)
>>> result = fit_density(x, basis=basis, lam=1e-2)
>>> grid = np.linspace(-5.0, 5.0, 2001)
>>> mass = float(np.sum(result(grid)) * (grid[1] - grid[0]))
>>> round(mass, 3)
1.0

DensityResult `dataclass` ¶

DensityResult(
    fd: FData,
    normaliser: float,
    lam: float,
    penalty_matrix: Array,
    criterion: float,
    log_likelihood: float,
    gradient_norm: float,
    n_iter: int,
    converged: bool,
)

Outcome of :func:fit_density.

The estimated density is p(x) = normaliser * exp(W(x)).

Attributes:

Name	Type	Description
`fd`	`FData`	The log-density function `W` up to the additive constant `log(normaliser)`. When `W` is only defined up to a constant (see the module notes) it is centred so that `∫ W = 0`.
`normaliser`	`float`	`C = 1 / ∫ exp W` over the basis domain.
`lam`	`float`	The smoothing parameter.
`penalty_matrix`	`array`	The roughness matrix `R`; the penalty is `lam * cᵀ R c`.
`criterion`	`float`	The minimised penalised negative log-likelihood `F(c)`.
`log_likelihood`	`float`	`Σ log p(x_i)` at the estimate.
`gradient_norm`	`float`	Euclidean norm of the gradient of `F` at the estimate.
`n_iter`	`int`	Number of Newton iterations taken.
`converged`	`bool`	Whether the step tolerance was met within `max_iter` iterations.

Examples:

>>> import numpy as np
>>> import fdatools as fdt
>>> from fdatools.density import fit_density
>>> x = np.array([0.1, 0.2, 0.25, 0.4, 0.45, 0.5, 0.55, 0.7, 0.8])
>>> res = fit_density(x, basis=fdt.BSpline(domain=(0.0, 1.0), n_basis=6), lam=1e-3)
>>> res.converged
True
>>> bool(np.all(res(np.array([0.0, 0.5, 1.0])) > 0))
True

call ¶

__call__(t: Any) -> Array

Evaluate the density p(t) = C exp W(t).

Parameters:

Name	Type	Description	Default
`t`	`array_like`	One-dimensional evaluation points inside the basis domain. A :class:`torch.Tensor` gives a tensor back.	required

Returns:

Type	Description
`array`	Density values, shape `(n_t,)`.

Examples:

>>> import numpy as np
>>> import fdatools as fdt
>>> from fdatools.density import fit_density
>>> x = np.array([0.2, 0.3, 0.5, 0.6, 0.65, 0.9])
>>> res = fit_density(x, basis=fdt.BSpline(domain=(0.0, 1.0), n_basis=5), lam=1.0)
>>> res(np.array([0.5])).shape
(1,)

log_density ¶

log_density(t: Any) -> Array

Evaluate the log-density log p(t) = W(t) + log C.

Parameters:

Name	Type	Description	Default
`t`	`array_like`	One-dimensional evaluation points inside the basis domain.	required

Returns:

Type	Description
`array`	Log-density values, shape `(n_t,)`.

Examples:

>>> import numpy as np
>>> import fdatools as fdt
>>> from fdatools.density import fit_density
>>> x = np.array([0.2, 0.3, 0.5, 0.6, 0.65, 0.9])
>>> res = fit_density(x, basis=fdt.BSpline(domain=(0.0, 1.0), n_basis=5), lam=1.0)
>>> t = np.array([0.25, 0.75])
>>> bool(np.allclose(np.exp(res.log_density(t)), res(t)))
True

IntensityResult `dataclass` ¶

IntensityResult(
    fd: FData,
    lam: float,
    penalty_matrix: Array,
    criterion: float,
    log_likelihood: float,
    expected_count: float,
    gradient_norm: float,
    n_iter: int,
    converged: bool,
)

Outcome of :func:fit_intensity.

The estimated intensity is μ(t) = exp(W(t)).

Attributes:

Name	Type	Description
`fd`	`FData`	The log-intensity function `W`.
`lam`	`float`	The smoothing parameter.
`penalty_matrix`	`array`	The roughness matrix `R`; the penalty is `lam * cᵀ R c`.
`criterion`	`float`	The minimised penalised negative log-likelihood `F(c)` (R's `Flist$f`).
`log_likelihood`	`float`	The Poisson-process log-likelihood `Σ W(x_i) - ∫ exp W`.
`expected_count`	`float`	`∫ exp W` over the domain, the expected number of events.
`gradient_norm`	`float`	Euclidean norm of the gradient of `F` at the estimate.
`n_iter`	`int`	Number of Newton iterations taken.
`converged`	`bool`	Whether the step tolerance was met within `max_iter` iterations.

Examples:

>>> import numpy as np
>>> import fdatools as fdt
>>> from fdatools.density import fit_intensity
>>> events = np.cumsum(np.full(40, 0.25))
>>> res = fit_intensity(events, basis=fdt.BSpline(domain=(0.0, 10.0), n_basis=5), lam=10.0)
>>> round(res.expected_count, 6)
40.0

call ¶

__call__(t: Any) -> Array

Evaluate the intensity μ(t) = exp W(t).

Parameters:

Name	Type	Description	Default
`t`	`array_like`	One-dimensional evaluation points inside the basis domain. A :class:`torch.Tensor` gives a tensor back.	required

Returns:

Type	Description
`array`	Intensity values (events per unit time), shape `(n_t,)`.

Examples:

>>> import numpy as np
>>> import fdatools as fdt
>>> from fdatools.density import fit_intensity
>>> events = np.linspace(0.5, 9.5, 19)
>>> basis = fdt.BSpline(domain=(0.0, 10.0), n_basis=4)
>>> res = fit_intensity(events, basis=basis, penalty=1, lam=1e6)
>>> np.round(res(np.array([2.0, 8.0])), 3).tolist()
[1.9, 1.9]

fit_density ¶

fit_density(
    x: Any,
    *,
    basis: Basis | None = None,
    domain: tuple[float, float] | None = None,
    lam: float = 0.0,
    penalty: int | LDO = 2,
    start: Any = None,
    tol: float = DEFAULT_TOL,
    max_iter: int = DEFAULT_MAX_ITER,
) -> DensityResult

Estimate a probability density by penalised maximum likelihood.

Replaces R fda's density.fd. The density is p(x) = C exp W(x) with W a basis expansion; W minimises -Σ W(x_i) + n log ∫ exp W + lam * ∫ (L W)².

Parameters:

Name	Type	Description	Default
`x`	`array_like`	One-dimensional sample (a :class:`torch.Tensor` is accepted and read as data; the fit itself is not differentiable).	required
`basis`	`Basis`	Basis for `W`; its domain is the support of the density. Default: cubic B-splines with `min(max(4, n // 10), 20)` functions on `domain`.	`None`
`domain`	`tuple of float`	Support used for the default basis; default `(min x, max x)`. Ignored when `basis` is given.	`None`
`lam`	`float`	Smoothing parameter `λ >= 0` (default `0`, no penalty).	`0.0`
`penalty`	`int or LDO`	Derivative order or linear differential operator `L` of the roughness penalty (default `2`). A heavy first-derivative penalty pulls the density towards the uniform one; a heavy `D³` penalty towards a (truncated) normal.	`2`
`start`	`FData or array_like`	Starting coefficients of `W`; default zero (the uniform density).	`None`
`tol`	`float`	Newton step tolerance: stop once a full Newton step moves no coefficient by more than `tol * (1 + max\|c\|)`.	`DEFAULT_TOL`
`max_iter`	`int`	Iteration cap.	`DEFAULT_MAX_ITER`

Returns:

Type	Description
`DensityResult`	`result(t)` evaluates the density, `result.fd` is `W` and `result.normaliser` is `C`.

Raises:

Type	Description
`ValueError`	If `x` is empty, not finite or not one-dimensional, lies outside the basis domain, or a setting is out of range.

Warns:

Type	Description
`RuntimeWarning`	If the Newton iteration has not converged after `max_iter` steps.

Notes

If the basis reproduces constants and L annihilates them, W is fixed only up to an additive constant; the returned W has ∫ W = 0. In that case the fit coincides with :func:fit_intensity on the same data, normalised: p = μ / n.

Examples:

>>> import numpy as np
>>> import fdatools as fdt
>>> from fdatools.density import fit_density
>>> rng = np.random.default_rng(1)
>>> x = rng.exponential(size=500)
>>> res = fit_density(x, domain=(0.0, 10.0), lam=1e-2)
>>> bool(res(np.array([0.1]))[0] > res(np.array([3.0]))[0])
True

fit_intensity ¶

fit_intensity(
    x: Any,
    *,
    basis: Basis | None = None,
    domain: tuple[float, float] | None = None,
    lam: float = 0.0,
    penalty: int | LDO = 2,
    start: Any = None,
    tol: float = DEFAULT_TOL,
    max_iter: int = DEFAULT_MAX_ITER,
) -> IntensityResult

Estimate the intensity of a Poisson process by penalised likelihood.

Replaces R fda's intensity.fd. The intensity is μ(t) = exp W(t); W minimises -Σ W(x_i) + ∫ exp W + lam * ∫ (L W)² over the basis domain.

Parameters:

Name	Type	Description	Default
`x`	`array_like`	One-dimensional event times (a :class:`torch.Tensor` is accepted and read as data). Their order does not matter.	required
`basis`	`Basis`	Basis for `W`; its domain is the observation window. Default: cubic B-splines with `min(max(4, n // 10), 20)` functions on `domain`.	`None`
`domain`	`tuple of float`	Observation window for the default basis; default `(0, max x)`, R's convention that observation starts at time 0. Ignored when `basis` is given.	`None`
`lam`	`float`	Smoothing parameter `λ >= 0` (default `0`).	`0.0`
`penalty`	`int or LDO`	Derivative order or operator `L` of the roughness penalty (default `2`). A heavy first-derivative penalty pulls the process towards a homogeneous one.	`2`
`start`	`FData or array_like`	Starting coefficients of `W`; default zero.	`None`
`tol`	`float`	Newton step tolerance: stop once a full Newton step moves no coefficient by more than `tol * (1 + max\|c\|)`.	`DEFAULT_TOL`
`max_iter`	`int`	Iteration cap.	`DEFAULT_MAX_ITER`

Returns:

Type	Description
`IntensityResult`	`result(t)` evaluates `μ(t)`; `result.fd` is `W`.

Raises:

Type	Description
`ValueError`	If `x` is empty, not finite or not one-dimensional, lies outside the basis domain, or a setting is out of range.

Warns:

Type	Description
`RuntimeWarning`	If the Newton iteration has not converged after `max_iter` steps.

Examples:

>>> import numpy as np
>>> import fdatools as fdt
>>> from fdatools.density import fit_intensity
>>> rng = np.random.default_rng(2)
>>> events = np.cumsum(rng.exponential(scale=0.5, size=200))
>>> res = fit_intensity(events, penalty=1, lam=10.0)
>>> round(res.expected_count, 6)
200.0

Density and intensity¶

DensityResult dataclass ¶

__call__ ¶

log_density ¶

IntensityResult dataclass ¶

__call__ ¶

fit_density ¶

fit_intensity ¶

DensityResult `dataclass` ¶

call ¶

IntensityResult `dataclass` ¶

call ¶