diff --git a/src/pysatl_core/distributions/support.py b/src/pysatl_core/distributions/support.py index 40f49bf..9aaf2e7 100644 --- a/src/pysatl_core/distributions/support.py +++ b/src/pysatl_core/distributions/support.py @@ -14,13 +14,25 @@ __copyright__ = "Copyright (c) 2025 PySATL project" __license__ = "SPDX-License-Identifier: MIT" +from collections.abc import Callable from dataclasses import dataclass from math import floor -from typing import TYPE_CHECKING, Protocol, cast, overload, runtime_checkable +from typing import ( + TYPE_CHECKING, + Protocol, + cast, + runtime_checkable, +) import numpy as np -from pysatl_core.types import BoolArray, Interval1D, Number, NumericArray +from pysatl_core.types import ( + BoolArray, + Interval1D, + IntervalND, + Number, + NumericArray, +) if TYPE_CHECKING: from collections.abc import Iterable, Iterator @@ -34,10 +46,7 @@ class Support(Protocol): Support defines the set of values where a distribution is defined. """ - @overload - def contains(self, x: Number) -> bool: ... - @overload - def contains(self, x: NumericArray) -> BoolArray: ... + def contains(self, x: NumericArray) -> bool | BoolArray: ... class ContinuousSupport(Interval1D, Support): @@ -49,6 +58,15 @@ class ContinuousSupport(Interval1D, Support): """ +class ContinuousNDSupport(IntervalND, Support): + """ + Support for continuous distributions represented as an array of intervals. + + This class inherits from IntervalND and implements the Support protocol + for continuous distributions defined on a list of intervals [left, right]. + """ + + @runtime_checkable class DiscreteSupport(Support, Protocol): """ @@ -128,12 +146,7 @@ def __init__(self, points: Iterable[Number], assume_sorted: bool = False) -> Non self._points = arr[unique_mask] - @overload - def contains(self, x: Number) -> bool: ... - @overload - def contains(self, x: NumericArray) -> BoolArray: ... - - def contains(self, x: Number | NumericArray) -> bool | BoolArray: + def contains(self, x: NumericArray) -> bool | BoolArray: """ Check if point(s) are in the support. @@ -162,10 +175,6 @@ def contains(self, x: Number | NumericArray) -> bool | BoolArray: return bool(result) return cast(BoolArray, result) - def __contains__(self, x: object) -> bool: - """Check if a point is in the support.""" - return bool(self.contains(cast(Number, x))) - def iter_points(self) -> Iterator[Number]: """Iterate through all points in the support.""" return iter(self._points) @@ -252,12 +261,7 @@ def __post_init__(self) -> None: if self.modulus <= 0: raise ValueError("modulus must be a positive integer.") - @overload - def contains(self, x: Number) -> bool: ... - @overload - def contains(self, x: NumericArray) -> BoolArray: ... - - def contains(self, x: Number | NumericArray) -> bool | BoolArray: + def contains(self, x: NumericArray) -> bool | BoolArray: """ Check if point(s) are in the integer lattice support. @@ -283,10 +287,6 @@ def contains(self, x: Number | NumericArray) -> bool | BoolArray: return bool(result) return cast(BoolArray, result) - def __contains__(self, x: object) -> bool: - """Check if a point is in the integer lattice support.""" - return bool(self.contains(cast(Number, x))) - def iter_points(self) -> Iterator[int]: """ Iterate through all points in the integer lattice support. @@ -430,10 +430,20 @@ def is_right_bounded(self) -> bool: __iter__ = iter_points +@dataclass(frozen=True, slots=True) +class PredicateSupport(Support): + predicate: Callable[[NumericArray], bool | BoolArray] + + def contains(self, x: NumericArray) -> bool | BoolArray: + return self.predicate(x) + + __all__ = [ # Base support protocol "Support", "ContinuousSupport", + "ContinuousNDSupport", + "PredicateSupport", # Discrete support protocol and implementations "DiscreteSupport", "ExplicitTableDiscreteSupport", diff --git a/src/pysatl_core/families/__init__.py b/src/pysatl_core/families/__init__.py index ed30528..502dac9 100644 --- a/src/pysatl_core/families/__init__.py +++ b/src/pysatl_core/families/__init__.py @@ -14,6 +14,12 @@ from .builtins import __all__ as _builtins_all from .configuration import configure_families_register from .distribution import ParametricFamilyDistribution +from .exponential_family import ( + # CanonicalContinuousExponentialClassFamily, + ContinuousExponentialClassFamily, + ExponentialConjugateHyperparameters, + ExponentialFamilyParametrization, +) from .parametric_family import ParametricFamily from .parametrizations import ( Parametrization, @@ -34,6 +40,9 @@ "configure_families_register", # builtins *_builtins_all, + "ContinuousExponentialClassFamily", + "ExponentialFamilyParametrization", + "ExponentialConjugateHyperparameters", ] del _builtins_all diff --git a/src/pysatl_core/families/builtins/continuous/__init__.py b/src/pysatl_core/families/builtins/continuous/__init__.py index d7f9492..66be601 100644 --- a/src/pysatl_core/families/builtins/continuous/__init__.py +++ b/src/pysatl_core/families/builtins/continuous/__init__.py @@ -11,10 +11,12 @@ from pysatl_core.families.builtins.continuous.exponential import configure_exponential_family from pysatl_core.families.builtins.continuous.normal import configure_normal_family +from pysatl_core.families.builtins.continuous.pareto import configure_pareto_family from pysatl_core.families.builtins.continuous.uniform import configure_uniform_family __all__ = [ "configure_normal_family", "configure_uniform_family", "configure_exponential_family", + "configure_pareto_family", ] diff --git a/src/pysatl_core/families/builtins/continuous/pareto.py b/src/pysatl_core/families/builtins/continuous/pareto.py new file mode 100644 index 0000000..5233979 --- /dev/null +++ b/src/pysatl_core/families/builtins/continuous/pareto.py @@ -0,0 +1,109 @@ +""" +Pareto distribution family implementation. + +This module provides a Pareto Type I distribution with fixed known minimum +``x_m = 1`` as a continuous exponential family. +""" + +from __future__ import annotations + +__author__ = "Vinogradov Ilya" +__copyright__ = "Copyright (c) 2025 PySATL project" +__license__ = "SPDX-License-Identifier: MIT" + + +import numpy as np + +from pysatl_core.distributions.support import ContinuousSupport +from pysatl_core.families.exponential_family import ( + ContinuousExponentialClassFamily, + ExponentialFamilyParametrization, +) +from pysatl_core.families.parametrizations import Parametrization, constraint, parametrization +from pysatl_core.families.registry import ParametricFamilyRegister +from pysatl_core.types import FamilyName, NumericArray, UnivariateContinuous + +_MINIMUM = 1.0 + + +def configure_pareto_family() -> None: + """ + Configure and register the Pareto distribution family with fixed minimum 1. + """ + + if ParametricFamilyRegister.contains(FamilyName.PARETO): + return + + PARETO_DOC = """ + Pareto Type I distribution with known minimum x_m = 1. + + The distribution is parameterized by the shape parameter ``alpha > 0``. + + Probability density function: + f(x) = alpha / x^(alpha + 1), x >= 1 + + In natural form this is written as + f(x | theta) = exp(theta log(x) + B(theta)), theta < -1 + where ``theta = -(alpha + 1)``. + + With the more common convention + f(x | theta) = exp(theta log(x) - A(theta)), + we have + A(theta) = -log(alpha * x_m^alpha) + and therefore + B(theta) = -A(theta) = log(alpha * x_m^alpha). + """ + + def _theta_to_alpha(theta: NumericArray) -> float: + theta_arr = np.atleast_1d(np.asarray(theta, dtype=float)) + return float(-theta_arr[0] - 1.0) + + def log_partition(theta: NumericArray) -> NumericArray: + alpha = _theta_to_alpha(theta) + return np.array([np.log(alpha) + alpha * np.log(_MINIMUM)]) + + def sufficient_statistics(x: NumericArray) -> NumericArray: + x_arr = np.asarray(x, dtype=float) + return np.array([np.log(x_arr).item()]) + + def normalization_constant(_: NumericArray) -> float: + return 1.0 + + def _support(_: Parametrization) -> ContinuousSupport: + return ContinuousSupport(left=_MINIMUM) + + pareto_family = ContinuousExponentialClassFamily( + name=FamilyName.PARETO, + log_partition=log_partition, + sufficient_statistics=sufficient_statistics, + normalization_constant=normalization_constant, + support=ContinuousSupport(left=_MINIMUM), + parameter_space=ContinuousSupport(right=-1.0, right_closed=False), + sufficient_statistics_values=ContinuousSupport(left=np.log(_MINIMUM)), + distr_type=UnivariateContinuous, + distr_parametrizations=["theta", "shape"], + support_by_parametrization=_support, + ) + pareto_family.__doc__ = PARETO_DOC + + @parametrization(family=pareto_family, name="shape") + class _Shape(Parametrization): + """ + Shape parametrization of Pareto distribution. + + Parameters + ---------- + alpha : float + Shape parameter of the distribution. + """ + + alpha: float + + @constraint(description="alpha > 0") + def check_alpha_positive(self) -> bool: + return self.alpha > 0 + + def transform_to_base_parametrization(self) -> ExponentialFamilyParametrization: + return ExponentialFamilyParametrization(theta=np.array([-(self.alpha + 1.0)])) + + ParametricFamilyRegister.register(pareto_family) diff --git a/src/pysatl_core/families/configuration.py b/src/pysatl_core/families/configuration.py index 7d0358a..5d5892b 100644 --- a/src/pysatl_core/families/configuration.py +++ b/src/pysatl_core/families/configuration.py @@ -26,6 +26,7 @@ from pysatl_core.families.builtins import ( configure_exponential_family, configure_normal_family, + configure_pareto_family, configure_uniform_family, ) from pysatl_core.families.registry import ParametricFamilyRegister @@ -46,6 +47,7 @@ def configure_families_register() -> ParametricFamilyRegister: The global registry of parametric families. """ configure_exponential_family() + configure_pareto_family() configure_uniform_family() configure_normal_family() return ParametricFamilyRegister() diff --git a/src/pysatl_core/families/exponential_family.py b/src/pysatl_core/families/exponential_family.py new file mode 100644 index 0000000..0169580 --- /dev/null +++ b/src/pysatl_core/families/exponential_family.py @@ -0,0 +1,540 @@ +""" +Exponential family distributions in continuous spaces. + +This module implements the continuous exponential family of probability distributions, +their conjugate priors, posterior inference, and posterior predictive distributions. +""" + +from __future__ import annotations + +__author__ = "Vinogradov Ilya" +__copyright__ = "Copyright (c) 2025 PySATL project" +__license__ = "SPDX-License-Identifier: MIT" + +from collections.abc import Callable, Iterable +from dataclasses import dataclass +from typing import TYPE_CHECKING, Any, cast + +import numpy as np +from scipy.differentiate import jacobian +from scipy.integrate import nquad +from scipy.linalg import det + +from pysatl_core.distributions.support import ( + ContinuousSupport, + PredicateSupport, +) +from pysatl_core.families.parametric_family import ParametricFamily +from pysatl_core.families.parametrizations import Parametrization, constraint, parametrization +from pysatl_core.types import ( + CharacteristicName, + DistributionType, + Number, + NumericArray, + ParametrizationName, + UnivariateContinuous, +) + +if TYPE_CHECKING: + from pysatl_core.distributions.support import Support + from pysatl_core.families.parametric_family import ( + CharacteristicsMap, + ParametricFamilyCharacteristic, + ) + + type SupportArg = Callable[[Parametrization], Support | None] | None + + +@dataclass +class ExponentialFamilyParametrization(Parametrization): + """ + Standard parametrization of an exponential family distribution. + + This parametrization uses the natural parameter vector ``theta``. In this + module the density is written as + + ``f(x | theta) = h(x) exp(theta^T T(x) + B(theta))``, + + where ``B(theta)`` is the log-normalizing term supplied as + ``log_partition``. + + Attributes + ---------- + theta : NumericArray + Natural parameter vector. + """ + + theta: NumericArray + + def transform_to_base_parametrization(self) -> ExponentialFamilyParametrization: + """Return the base parametrization (identity transform for canonical form).""" + return self + + +@dataclass +class ExponentialConjugateHyperparameters(Parametrization): + """ + Hyperparameters for the conjugate prior of an exponential family. + + For this module's sign convention, the conjugate prior over ``theta`` is + proportional to + + ``exp(nu_0^T theta + n_0 B(theta))``, + + where ``B(theta)`` is the base family's log-normalizing term. + + Attributes + ---------- + effective_suff_stat_value : NumericArray + Pseudo-sufficient statistic value ``nu_0``. + effective_sample_size : Number + Pseudo-sample size ``n_0``. + """ + + effective_suff_stat_value: NumericArray + effective_sample_size: Number + + def transform_to_base_parametrization(self) -> ExponentialFamilyParametrization: + """ + Convert hyperparameters to a canonical parametrization. + + The resulting parameter vector is [ν₀, n₀] concatenated. + """ + return ExponentialFamilyParametrization( + np.append(self.effective_suff_stat_value, self.effective_sample_size) + ) + + +class ContinuousExponentialClassFamily(ParametricFamily): + """ + Representation of a continuous exponential family distribution. + + The density is given by + + ``f(x | theta) = h(x) exp(theta^T T(x) + B(theta))``. + + Here ``theta`` is the natural parameter, ``T(x)`` is the sufficient + statistic, ``h(x)`` is the base measure supplied as + ``normalization_constant``, and ``B(theta)`` is the log-normalizing term + supplied as ``log_partition``. With the usual convention + ``h(x) exp(theta^T T(x) - A(theta))``, this means ``B(theta) = -A(theta)``. + + The family provides canonical parametrization by ``theta``, conjugate prior + construction, posterior hyperparameter updates, posterior predictive + densities, and monotone differentiable transformations. + """ + + def __init__( + self, + *, + log_partition: Callable[[NumericArray], NumericArray], + sufficient_statistics: Callable[[NumericArray], NumericArray], + normalization_constant: Callable[[NumericArray], Number], + support: Support, + parameter_space: Support, + sufficient_statistics_values: Support, + name: str, + distr_type: DistributionType | Callable[[Parametrization], DistributionType], + distr_parametrizations: list[ParametrizationName], + distr_characteristics: CharacteristicsMap | None = None, + support_by_parametrization: SupportArg = None, + base_score: Callable[[Parametrization, NumericArray], NumericArray] | None = None, + ): + """ + Initialize a continuous exponential family distribution. + + Parameters + ---------- + log_partition : Callable[[NumericArray], NumericArray] + Function ``B(theta)`` used in + ``f(x | theta) = h(x) exp(theta^T T(x) + B(theta))``. + sufficient_statistics : Callable[[NumericArray], NumericArray] + Function ``T(x)`` returning the sufficient statistic. + normalization_constant : Callable[[NumericArray], Number] + Function ``h(x)`` returning the base measure. + support : Support + Support of the observation variable ``x``. + parameter_space : Support + Support of the natural parameter ``theta``. + sufficient_statistics_values : Support + Support of possible sufficient-statistic values. + name : str + Family name. + distr_type : DistributionType or Callable[[Parametrization], DistributionType] + Distribution type, or a resolver from base parametrization to type. + distr_parametrizations : list[ParametrizationName] + Parametrization names supported by this family. + distr_characteristics : CharacteristicsMap, optional + Additional analytical characteristics to register. + support_by_parametrization : Callable[[Parametrization], Support | None], optional + Resolver for the distribution support of a concrete parametrization. + base_score : Callable[[Parametrization, NumericArray], NumericArray], optional + Score function in the base parametrization. + """ + self._sufficient = sufficient_statistics + self._log_partition = log_partition + self._normalization = normalization_constant + + self._support = support + self._parameter_space = parameter_space + self._sufficient_statistics_values = sufficient_statistics_values + + family_characteristics: CharacteristicsMap = { + CharacteristicName.PDF: self.density, + CharacteristicName.MEAN: self._mean, + CharacteristicName.VAR: self._var, + } + distr_characteristics = dict(distr_characteristics or {}) + merged_characteristics = dict(family_characteristics) + merged_characteristics.update(distr_characteristics) + + ParametricFamily.__init__( + self, + name=name, + distr_type=distr_type, + distr_parametrizations=distr_parametrizations, + distr_characteristics=merged_characteristics, + support_by_parametrization=support_by_parametrization, + base_score=base_score, + ) + + family = self + + @parametrization(family=family, name="theta") + class ThetaParametrization(ExponentialFamilyParametrization): + @constraint(description="theta belongs to parameter_space") + def check_theta_in_parameter_space(self) -> bool: + theta = np.atleast_1d(np.asarray(self.theta, dtype=float)) + return family._parameter_space_contains(theta) + + def _parameter_space_contains(self, theta: NumericArray) -> bool: + is_contains = self._parameter_space.contains(theta) + if isinstance(is_contains, np.ndarray): + return bool(np.all(is_contains)) + return bool(is_contains) + + @property + def log_density(self) -> ParametricFamilyCharacteristic[NumericArray, Number]: + """ + Return the log-density characteristic. + + The returned callable evaluates + + ``log h(x) + theta^T T(x) + B(theta)``. + + Points outside the observation support return ``-np.inf``. + + Returns + ------- + ParametricFamilyCharacteristic[NumericArray, Number] + Callable accepting a parametrization and observation. + """ + + def log_density_func(parametrization: Parametrization, x: NumericArray) -> Number: + parametrization = cast(ExponentialFamilyParametrization, parametrization) + parametrization = parametrization.transform_to_base_parametrization() + if not self._support.contains(np.array([x])): + return -np.inf + + theta = parametrization.theta + sufficient = self._sufficient(x) + dot = np.dot(theta, sufficient) + if hasattr(dot, "__len__"): + dot = dot[0] + + result = np.log(self._normalization(x)) + dot + self._log_partition(theta) + return cast(np.floating, result.item()) + + return log_density_func + + @property + def density(self) -> ParametricFamilyCharacteristic[NumericArray, Number]: + """ + Return the density characteristic. + + Returns + ------- + ParametricFamilyCharacteristic[NumericArray, Number] + Callable evaluating ``exp(log_density(parametrization, x))``. + """ + log_density = cast(Callable[[Parametrization, NumericArray], Number], self.log_density) + + def density_func(parametrization: Parametrization, x: NumericArray) -> Number: + return cast(Number, np.exp(log_density(parametrization, x))) + + return density_func + + @property + def conjugate_prior_family(self) -> ContinuousExponentialClassFamily: + """ + Build the conjugate prior family for this exponential family. + + The conjugate prior is an exponential family over ``theta`` with + sufficient statistic ``[theta, B(theta)]`` and base measure ``1``. + Its natural parameter is ``[nu, n]``, matching + :class:`ExponentialConjugateHyperparameters`. + + Returns + ------- + ContinuousExponentialClassFamily + Conjugate prior family for the natural parameter. + """ + + def conjugate_sufficient( + theta: NumericArray, + ) -> NumericArray: + if not hasattr(theta, "__len__"): + theta = np.array([theta]) + + if not self._parameter_space.contains(theta): + return np.full(len(theta) + 1, float("-inf")) + return np.append(theta, self._log_partition(theta)) + + def conjugate_log_partition( + parametrization: NumericArray, + ) -> NumericArray: + def pdf(theta: NumericArray) -> Number: + if not hasattr(theta, "__len__"): + theta = np.array([theta]) + return cast( + np.floating, + np.exp( + np.dot( + conjugate_sufficient(theta), + parametrization, + ) + ).item(), + ) + + def integrand(x: float) -> float: + theta = np.asarray([x], dtype=float) + if not self._parameter_space.contains(theta): + return 0.0 + return float(pdf(theta)) + + all_value = nquad(integrand, [(float("-inf"), float("+inf"))])[0] + return np.array([cast(np.float64, -np.log(all_value))]) + + def conjugate_sufficient_accepts( + theta: NumericArray, + ) -> bool: + xi = theta[:-1] + nu = theta[-1] + + return bool(self._sufficient_statistics_values.contains(xi)) and bool( + ContinuousSupport(0, np.inf).contains(np.array([nu])) + ) + + return ContinuousExponentialClassFamily( + log_partition=conjugate_log_partition, + sufficient_statistics=conjugate_sufficient, + normalization_constant=lambda _: 1, + support=self._parameter_space, + sufficient_statistics_values=self._parameter_space, + parameter_space=PredicateSupport(predicate=conjugate_sufficient_accepts), + name=self.name, + distr_type=self._distr_type, + distr_parametrizations=self.parametrization_names, + support_by_parametrization=self.support_resolver, + ) + + def transform( + self, + inverse_transform_function: Callable[[NumericArray], NumericArray], + ) -> ContinuousExponentialClassFamily: + """ + Transform the random variable by a monotonic, differentiable function. + + ``transform_function`` is interpreted as the inverse map ``x = g(y)``. + The transformed density is + + ``f_Y(y | theta) = h(g(y)) exp(theta^T T(g(y)) + B(theta)) |J_g(y)|``. + + Parameters + ---------- + transform_function : Callable[[NumericArray], NumericArray] + Invertible differentiable function ``g`` mapping transformed + observations back to the original observation support. + + Returns + ------- + ContinuousExponentialClassFamily + Family for the transformed random variable. + """ + + def calculate_jacobian(x: NumericArray) -> NumericArray: + if not isinstance(x, Iterable): + x = np.array([x], dtype=float) + else: + x = np.atleast_1d(np.asarray(x, dtype=float)) + + return np.abs(det(jacobian(inverse_transform_function, x).df)) + + def new_support(x: NumericArray) -> bool: + return bool(self._support.contains(np.asarray(inverse_transform_function(x)))) + + def new_sufficient(x: NumericArray) -> NumericArray: + return self._sufficient(inverse_transform_function(x)) + + def new_normalization(x: NumericArray) -> Number: + return cast( + np.float64, + self._normalization(inverse_transform_function(x)) * calculate_jacobian(x), + ) + + return ContinuousExponentialClassFamily( + log_partition=self._log_partition, + sufficient_statistics=new_sufficient, + normalization_constant=new_normalization, + support=PredicateSupport(predicate=new_support), + parameter_space=self._parameter_space, + sufficient_statistics_values=self._sufficient_statistics_values, + name=f"Transformed{self._name}", + distr_type=self._distr_type, + distr_parametrizations=self.parametrization_names, + support_by_parametrization=self.support_resolver, + ) + + @property + def _mean(self) -> ParametricFamilyCharacteristic[Any, Any]: + """Compute the mean E[X] by numerical integration over the density.""" + + def mean_func(parametrization: Parametrization) -> Any: + parametrization = cast(ExponentialFamilyParametrization, parametrization) + density = cast(Callable[[Parametrization, NumericArray], Number], self.density) + return nquad( + lambda x: ( + np.dot(x, density(parametrization, x)) + if self._support.contains(np.array([x])) + else 0 + ), + [(float("-inf"), float("inf"))], + )[0] + + return mean_func + + @property + def _second_moment(self) -> ParametricFamilyCharacteristic[Any, Any]: + """Compute the second moment E[X²] by numerical integration.""" + + def func(parametrization: Parametrization) -> Any: + parametrization = cast(ExponentialFamilyParametrization, parametrization) + density = cast(Callable[[Parametrization, NumericArray], Number], self.density) + return nquad( + lambda x: ( + x**2 * density(parametrization, x) + if self._support.contains(np.array([x])) + else 0 + ), + [(float("-inf"), float("inf"))], + )[0] + + return func + + @property + def _var(self) -> ParametricFamilyCharacteristic[Any, Any]: + """Compute the variance Var[X] = E[X²] - (E[X])².""" + + def func(parametrization: Parametrization) -> Any: + parametrization = cast(ExponentialFamilyParametrization, parametrization) + second_moment = cast(Callable[[Parametrization], Any], self._second_moment) + mean = cast(Callable[[Parametrization], Any], self._mean) + return second_moment(parametrization) - mean(parametrization) ** 2 + + return func + + def posterior_hyperparameters( + self, + parametrization: ExponentialConjugateHyperparameters, + sample: list[Any] | Any, + ) -> ExponentialConjugateHyperparameters: + """ + Update the conjugate prior hyperparameters given observed data. + + For prior hyperparameters ``(nu_0, n_0)`` and observations ``x_i``, + the posterior hyperparameters are + + ``nu = nu_0 + sum_i T(x_i)`` and ``n = n_0 + N``. + + Parameters + ---------- + parametrization : ExponentialConjugateHyperparameters + Current conjugate hyperparameters. + sample : list[Any] or Any + Observations used for the update. A non-string iterable is treated + as a sample; any other value is treated as one observation. + + Returns + ------- + ExponentialConjugateHyperparameters + Updated hyperparameters after incorporating ``sample``. + """ + if hasattr(sample, "__iter__") and not isinstance(sample, str): + posterior_effective_suff_stat_value = np.array( + parametrization.effective_suff_stat_value + ) + np.sum( + [self._sufficient(x) for x in sample], + axis=0, + ) + posterior_effective_sample_size = parametrization.effective_sample_size + len(sample) + else: + posterior_effective_suff_stat_value = np.array( + parametrization.effective_suff_stat_value, + ) + np.asarray(self._sufficient(sample)) # type: ignore[arg-type] + posterior_effective_sample_size = parametrization.effective_sample_size + 1 + + return ExponentialConjugateHyperparameters( + effective_suff_stat_value=posterior_effective_suff_stat_value, + effective_sample_size=posterior_effective_sample_size, + ) + + @property + def posterior_predictive(self) -> ParametricFamily: + """ + Construct the posterior predictive distribution. + + For conjugate hyperparameters ``(nu, n)``, the predictive density of a + new observation ``x`` is + + ``p(x | nu, n) = h(x) Z_c(nu, n) / Z_c(nu + T(x), n + 1)``, + + where ``Z_c`` is the conjugate-prior normalizing integral. The density + is zero outside the original observation support. + + Returns + ------- + ParametricFamily + Family with ``ExponentialConjugateHyperparameters`` parametrization + and a ``pdf`` characteristic for the posterior predictive density. + """ + + def conjugate_log_partition( + parametrization: ExponentialConjugateHyperparameters, + ) -> Number: + conjugate_value = self.conjugate_prior_family._log_partition( + parametrization.transform_to_base_parametrization().theta + ) + return cast(Number, np.exp(conjugate_value).item()) + + def posterior_density(parametrization: Parametrization, x: NumericArray) -> Number: + parametrization = cast(ExponentialConjugateHyperparameters, parametrization) + if not self._support.contains(np.asarray(x)): + return cast(np.float32, 0.0) + return cast( + np.float32, + self._normalization(x) + * conjugate_log_partition(parametrization) + / conjugate_log_partition( + self.posterior_hyperparameters(parametrization=parametrization, sample=[x]) + ), + ) + + family = ParametricFamily( + name=f"PosteriorPredictive{self.name}", + distr_type=UnivariateContinuous, + distr_characteristics={CharacteristicName.PDF: posterior_density}, + distr_parametrizations=["posterior"], + support_by_parametrization=lambda _: self._support, + ) + parametrization(family=family, name="posterior")(ExponentialConjugateHyperparameters) + return family diff --git a/src/pysatl_core/types.py b/src/pysatl_core/types.py index 25aa63e..de4a735 100644 --- a/src/pysatl_core/types.py +++ b/src/pysatl_core/types.py @@ -12,7 +12,7 @@ from dataclasses import dataclass from enum import Enum, StrEnum, auto from math import inf -from typing import TYPE_CHECKING, Any, cast, overload +from typing import TYPE_CHECKING, Any if TYPE_CHECKING: from pysatl_core.distributions.computations.computation import ( @@ -179,13 +179,7 @@ def __post_init__(self) -> None: if self.right == inf and self.right_closed: object.__setattr__(self, "right_closed", False) - @overload - def contains(self, x: Number) -> bool: ... - - @overload - def contains(self, x: NumericArray) -> BoolArray: ... - - def contains(self, x: Number | NumericArray) -> bool | BoolArray: + def contains(self, x: NumericArray) -> bool | BoolArray: """ Check if point(s) are contained in the interval. @@ -210,10 +204,6 @@ def contains(self, x: Number | NumericArray) -> bool | BoolArray: return result - def __contains__(self, x: object) -> bool: - """Check if a single point is in the interval.""" - return bool(self.contains(cast(Number, x))) - @property def is_empty(self) -> bool: """Check if the interval is empty.""" @@ -250,6 +240,25 @@ def shape(self) -> ContinuousSupportShape1D: type Method[In, Out] = AnalyticalComputation[In, Out] | FittedComputationMethod[In, Out] """Type alias for a distribution computation method (analytical or fitted).""" + +@dataclass(frozen=True, slots=True) +class IntervalND: + intervals: list[Interval1D] + + def contains(self, x: NumericArray) -> bool | BoolArray: + def contains_for_point(point: NumericArray) -> bool: + assert len(point) == len(self.intervals) + return all( + bool(interval.contains(np.asarray(x_coordinate))) + for interval, x_coordinate in zip(self.intervals, point, strict=True) + ) + + if len(x.shape) == 1: + return contains_for_point(x) + + return np.array([contains_for_point(point) for point in x]) + + type GenericCharacteristicName = str """Type alias for characteristic names (e.g., 'pdf', 'cdf').""" @@ -404,6 +413,7 @@ class FamilyName(StrEnum): NORMAL = "Normal" CONTINUOUS_UNIFORM = "ContinuousUniform" EXPONENTIAL = "Exponential" + PARETO = "Pareto" # ============================================================================ @@ -465,6 +475,7 @@ class FamilyName(StrEnum): "TransformationMethodSpecsMap", "DistributionType", "Interval1D", + "IntervalND", "ContinuousSupportShape1D", "BoolArray", "NumPyNumber", diff --git a/tests/unit/distributions/test_support.py b/tests/unit/distributions/test_support.py index 11d70bf..c6d5ceb 100644 --- a/tests/unit/distributions/test_support.py +++ b/tests/unit/distributions/test_support.py @@ -40,16 +40,14 @@ class TestContinuousSupport: ], ) def test_continuous_support_contains_scalar(self, point, expected_result): - assert (point in self.support_example) is expected_result - assert self.support_example.contains(point) is expected_result + assert self.support_example.contains(np.asarray(point)) is expected_result @pytest.mark.parametrize("infinity", [-inf, inf]) def test_continuous_support_doesnt_contain_inf(self, infinity): # inf isn't considered as a number # but as a limit so support doesn't contain it even if it's a real line support = ContinuousSupport() - assert infinity not in support - assert support.contains(infinity) is False + assert support.contains(np.asarray(infinity)) is False @pytest.mark.parametrize( "points,expected_result", @@ -59,7 +57,6 @@ def test_continuous_support_doesnt_contain_inf(self, infinity): ], ) def test_continuous_support_contains_array(self, points, expected_result): - # np.array doesn't have `in`(__contains__) syntax result = self.support_example.contains(points) assert isinstance(result, np.ndarray) assert result.tolist() == expected_result @@ -104,8 +101,7 @@ def test_table_is_sorted_and_deduplicated(self): ], ) def test_contains_scalar(self, point, expected_result): - assert (point in self.support_example) is expected_result - assert self.support_example.contains(point) is expected_result + assert self.support_example.contains(np.asarray(point)) is expected_result @pytest.mark.parametrize( "points, expected_result", @@ -113,7 +109,7 @@ def test_contains_scalar(self, point, expected_result): (np.array([0, 1, 2, 3, 4, 5]), [False, True, True, True, False, True]), (np.array([]), []), (np.array([1.0, 1.5, 2.0, 2.5]), [True, False, True, False]), - ([0, 1, 2, 3, 4, 5], [False, True, True, True, False, True]), + (np.array([0, 1, 2, 3, 4, 5]), [False, True, True, True, False, True]), ], ) def test_contains_array(self, points, expected_result): @@ -196,8 +192,7 @@ def test_invalid_modulus_raises(self): ) def test_contains_scalar(self, support_name, point, expected_result): support = self.support_examples[support_name] - assert (point in support) is expected_result - assert support.contains(point) is expected_result + assert support.contains(np.asarray(point)) is expected_result @pytest.mark.parametrize( "support_name, points, expected_result", diff --git a/tests/unit/families/builtins/continuous/test_exponential.py b/tests/unit/families/builtins/continuous/test_exponential.py index 6363c47..e92ce65 100644 --- a/tests/unit/families/builtins/continuous/test_exponential.py +++ b/tests/unit/families/builtins/continuous/test_exponential.py @@ -206,10 +206,10 @@ def test_exponential_support(self): assert not dist.support.right_closed # Test containment - assert dist.support.contains(0.0) is True - assert dist.support.contains(1.0) is True - assert dist.support.contains(-0.1) is False - assert dist.support.contains(float("inf")) is False + assert dist.support.contains(np.asarray(0.0)) is True + assert dist.support.contains(np.asarray(1.0)) is True + assert dist.support.contains(np.asarray(-0.1)) is False + assert dist.support.contains(np.asarray(float("inf"))) is False # Test array test_points = np.array([-0.1, 0.0, 1.0, 10.0]) diff --git a/tests/unit/families/builtins/continuous/test_normal.py b/tests/unit/families/builtins/continuous/test_normal.py index ca59349..0f47a42 100644 --- a/tests/unit/families/builtins/continuous/test_normal.py +++ b/tests/unit/families/builtins/continuous/test_normal.py @@ -226,9 +226,9 @@ def test_normal_support(self): assert not dist.support.left_closed assert not dist.support.right_closed - assert dist.support.contains(0) is True - assert dist.support.contains(float("inf")) is False - assert dist.support.contains(float("-inf")) is False + assert dist.support.contains(np.asarray(0)) is True + assert dist.support.contains(np.asarray(float("inf"))) is False + assert dist.support.contains(np.asarray(float("-inf"))) is False test_points = np.array([-500, 0, 5]) results = dist.support.contains(test_points) diff --git a/tests/unit/families/builtins/continuous/test_pareto.py b/tests/unit/families/builtins/continuous/test_pareto.py new file mode 100644 index 0000000..d362dcc --- /dev/null +++ b/tests/unit/families/builtins/continuous/test_pareto.py @@ -0,0 +1,168 @@ +""" +Tests for Pareto distribution family. +""" + +__author__ = "Vinogradov Ilya" +__copyright__ = "Copyright (c) 2025 PySATL project" +__license__ = "SPDX-License-Identifier: MIT" + +from typing import Any, cast + +import numpy as np +import pytest +from numpy.testing import assert_allclose +from scipy.stats import pareto + +from pysatl_core.distributions.support import ContinuousSupport +from pysatl_core.families.configuration import configure_families_register +from pysatl_core.families.exponential_family import ( + ContinuousExponentialClassFamily, + ExponentialConjugateHyperparameters, +) +from pysatl_core.types import ( + CharacteristicName, + ContinuousSupportShape1D, + FamilyName, + UnivariateContinuous, +) + +from .base import BaseDistributionTest + + +class TestParetoFamily(BaseDistributionTest): + """Test suite for Pareto distribution family.""" + + def setup_method(self) -> None: + registry = configure_families_register() + self.pareto_family = cast(ContinuousExponentialClassFamily, registry.get(FamilyName.PARETO)) + self.pareto_dist_example = self.pareto_family(parametrization_name="shape", alpha=3.0) + + def test_family_properties(self) -> None: + assert self.pareto_family.name == FamilyName.PARETO + assert isinstance(self.pareto_family, ContinuousExponentialClassFamily) + assert self.pareto_family.distribution( + parametrization_name="shape", alpha=3.0 + ).family_name == (FamilyName.PARETO) + assert set(self.pareto_family.parametrization_names) == {"theta", "shape"} + assert self.pareto_family.base_parametrization_name == "theta" + + def test_shape_parametrization_creation(self) -> None: + dist = self.pareto_family(parametrization_name="shape", alpha=3.0) + + assert dist.distribution_type == UnivariateContinuous + assert dist.parameters == {"alpha": 3.0} + assert dist.parametrization_name == "shape" + + def test_parametrization_constraints(self) -> None: + with pytest.raises(ValueError, match="alpha > 0"): + self.pareto_family(parametrization_name="shape", alpha=0.0) + + def test_parametrization_conversion_to_theta(self) -> None: + shape_params = cast(Any, self.pareto_family.get_parametrization("shape")) + base_params = cast( + Any, + self.pareto_family.to_base(shape_params(alpha=3.0)), + ) + + assert_allclose(base_params.theta, np.array([-4.0])) + + def test_analytical_computations_availability(self) -> None: + comp = self.pareto_dist_example.analytical_computations + + expected_chars = { + CharacteristicName.PDF, + CharacteristicName.MEAN, + CharacteristicName.VAR, + } + assert set(comp.keys()) == expected_chars + + @pytest.mark.parametrize( + "char_name, test_data, scipy_func", + [ + (CharacteristicName.PDF, [0.5, 1.0, 1.5, 2.0, 4.0], pareto.pdf), + ], + ) + def test_main_characteristics_against_scipy( + self, + char_name: CharacteristicName, + test_data: list[float], + scipy_func: Any, + ) -> None: + char_func = self.pareto_dist_example.query_method(char_name) + input_array = np.array(test_data) + if char_name == CharacteristicName.PDF: + actual = np.asarray([char_func(x) for x in input_array], dtype=float) + else: + actual = char_func(input_array) + expected = scipy_func(input_array, 3.0, scale=1.0) + + assert actual.shape == input_array.shape + self.assert_arrays_almost_equal(actual, expected, precision=1e-8) + + def test_log_pdf_matches_log_of_pdf(self) -> None: + pdf = self.pareto_dist_example.query_method(CharacteristicName.PDF) + x_values = np.array([1.0, 1.5, 2.0, 4.0]) + + actual = np.log(np.asarray([pdf(x) for x in x_values], dtype=float)) + expected = pareto.logpdf(x_values, 3.0, scale=1.0) + + self.assert_arrays_almost_equal(actual, expected, precision=1e-8) + + def test_moments(self) -> None: + mean_func = self.pareto_dist_example.query_method(CharacteristicName.MEAN) + var_func = self.pareto_dist_example.query_method(CharacteristicName.VAR) + + scipy_pareto = pareto(b=3.0, scale=1.0) + + assert mean_func() == pytest.approx(float(scipy_pareto.mean()), rel=1e-6) + assert var_func() == pytest.approx(float(scipy_pareto.var()), rel=1e-6) + + def test_support(self) -> None: + dist = self.pareto_dist_example + + assert dist.support is not None + assert isinstance(dist.support, ContinuousSupport) + assert dist.support.left == 1.0 + assert dist.support.right == float("inf") + assert dist.support.left_closed + assert not dist.support.right_closed + assert dist.support.shape == ContinuousSupportShape1D.RAY_RIGHT + + assert dist.support.contains(np.asarray(1.0)) is True + assert dist.support.contains(np.asarray(2.0)) is True + assert dist.support.contains(np.asarray(0.99)) is False + + def test_posterior_hyperparameters(self) -> None: + predictive_family = self.pareto_family.posterior_predictive + posterior_params_cls = cast(Any, predictive_family.get_parametrization("posterior")) + prior = cast( + ExponentialConjugateHyperparameters, + posterior_params_cls( + effective_suff_stat_value=np.array([2.0]), + effective_sample_size=3.0, + ), + ) + posterior = self.pareto_family.posterior_hyperparameters(prior, sample=[2.0, 4.0]) + + assert_allclose( + posterior.effective_suff_stat_value, + np.array([2.0 + np.log(2.0) + np.log(4.0)]), + ) + assert posterior.effective_sample_size == 5.0 + + def test_posterior_predictive_pdf(self) -> None: + xi = 2.0 + nu = 3.0 + predictive = self.pareto_family.posterior_predictive( + parametrization_name="posterior", + effective_suff_stat_value=np.array([xi]), + effective_sample_size=nu, + ) + pdf = predictive.query_method(CharacteristicName.PDF) + x_values = np.array([1.0, 1.5, 2.0, 4.0]) + + actual = np.asarray([pdf(x) for x in x_values], dtype=float) + expected = ( + (nu + 1.0) * xi ** (nu + 1.0) / (x_values * (xi + np.log(x_values)) ** (nu + 2.0)) + ) + self.assert_arrays_almost_equal(actual, expected, precision=1e-6) diff --git a/tests/unit/families/builtins/continuous/test_uniform.py b/tests/unit/families/builtins/continuous/test_uniform.py index 8d63f7f..d0cffc3 100644 --- a/tests/unit/families/builtins/continuous/test_uniform.py +++ b/tests/unit/families/builtins/continuous/test_uniform.py @@ -230,11 +230,11 @@ def test_uniform_support(self): assert dist.support.right_closed # Test containment - assert dist.support.contains(2.0) is True - assert dist.support.contains(5.0) is True - assert dist.support.contains(3.5) is True - assert dist.support.contains(1.9) is False - assert dist.support.contains(5.1) is False + assert dist.support.contains(np.asarray(2.0)) is True + assert dist.support.contains(np.asarray(5.0)) is True + assert dist.support.contains(np.asarray(3.5)) is True + assert dist.support.contains(np.asarray(1.9)) is False + assert dist.support.contains(np.asarray(5.1)) is False # Test array test_points = np.array([1.9, 2.0, 3.5, 5.0, 5.1]) diff --git a/tests/unit/families/test_configuration.py b/tests/unit/families/test_configuration.py index 0a7f2d8..546582a 100644 --- a/tests/unit/families/test_configuration.py +++ b/tests/unit/families/test_configuration.py @@ -40,6 +40,7 @@ def test_families_registered(self): expected_families = { FamilyName.NORMAL, FamilyName.CONTINUOUS_UNIFORM, + FamilyName.PARETO, } registered_families = set(self.registry._registered_families.keys()) @@ -78,4 +79,5 @@ def test_registry_list_registered_families(self): assert isinstance(families_list, list) assert FamilyName.NORMAL in families_list assert FamilyName.CONTINUOUS_UNIFORM in families_list + assert FamilyName.PARETO in families_list assert "NonExistentFamily" not in families_list diff --git a/tests/unit/families/test_exponential_family.py b/tests/unit/families/test_exponential_family.py new file mode 100644 index 0000000..b1e4ef1 --- /dev/null +++ b/tests/unit/families/test_exponential_family.py @@ -0,0 +1,333 @@ +from collections.abc import Callable + +__author__ = "Vinogradov Ilya" +__copyright__ = "Copyright (c) 2025 PySATL project" +__license__ = "SPDX-License-Identifier: MIT" + +import itertools +from typing import cast + +import numpy as np +import pytest +import scipy +from numpy.testing import assert_allclose + +from pysatl_core.distributions.support import ( + ContinuousNDSupport, + PredicateSupport, +) +from pysatl_core.families import ( + ContinuousExponentialClassFamily, + ExponentialConjugateHyperparameters, + ExponentialFamilyParametrization, +) +from pysatl_core.families.registry import ParametricFamilyRegister +from pysatl_core.types import ( + CharacteristicName, + Interval1D, + Number, + NumericArray, + UnivariateContinuous, +) + + +def gamma_pdf(alpha: float, beta: float, x: float) -> float: + return scipy.stats.gamma(a=alpha, scale=1 / beta).pdf(x).item() # type: ignore[attr-defined] + + +def lomax_pdf(shape: float, scale: float, x: float) -> float: + return scipy.stats.lomax(c=shape, scale=scale).pdf(x).item() # type: ignore[attr-defined] + + +def exponential_log_partition(parametrization): + return np.log(-parametrization) + + +def _make_exponential_family( + normalization_constant: Callable[[NumericArray], Number] = lambda _: 1, +) -> ContinuousExponentialClassFamily: + support_neg = PredicateSupport( + predicate=lambda x: bool( + ContinuousNDSupport( + intervals=[Interval1D(-np.inf, 0, left_closed=False, right_closed=False)] + ).contains(np.array([x])) + ) + ) + support_pos = PredicateSupport( + predicate=lambda x: bool( + ContinuousNDSupport( + intervals=[Interval1D(0, np.inf, left_closed=False, right_closed=False)] + ).contains(np.array([x])) + ) + ) + return ContinuousExponentialClassFamily( + name="ExponentialFamily", + log_partition=exponential_log_partition, + sufficient_statistics=lambda x: x, + normalization_constant=normalization_constant, + parameter_space=support_neg, + sufficient_statistics_values=support_pos, + support=support_pos, + distr_type=UnivariateContinuous, + distr_parametrizations=["theta"], + ) + + +@pytest.fixture(scope="function") +def exponential_family() -> ContinuousExponentialClassFamily: + return _make_exponential_family() + + +@pytest.fixture(scope="function") +def conjugate_for_exponential() -> ContinuousExponentialClassFamily: + def transform_function(x: NumericArray) -> NumericArray: + return -x + + fam = _make_exponential_family() + conjugate_fam = fam.conjugate_prior_family.transform(transform_function) + ParametricFamilyRegister().register(conjugate_fam) + return cast( + ContinuousExponentialClassFamily, + ParametricFamilyRegister().get("TransformedExponentialFamily"), + ) + + +def test_log_density_matches_exponential_form( + exponential_family: ContinuousExponentialClassFamily, +) -> None: + params = ExponentialFamilyParametrization(theta=np.array([-2.0])) + log_density_func = cast( + Callable[[ExponentialFamilyParametrization, NumericArray], Number], + exponential_family.log_density, + ) + density_func = cast( + Callable[[ExponentialFamilyParametrization, NumericArray], Number], + exponential_family.density, + ) + + log_density = log_density_func(params, np.asarray(0.5)) + density = density_func(params, np.asarray(0.5)) + + assert log_density == pytest.approx(np.log(2.0) - 1.0) + assert density == pytest.approx(2.0 * np.exp(-1.0)) + + +def test_constructor_merges_custom_characteristics() -> None: + support_neg = PredicateSupport( + predicate=lambda x: bool( + ContinuousNDSupport(intervals=[Interval1D(-np.inf, 0)]).contains(np.array([x])) + ) + ) + support_pos = PredicateSupport( + predicate=lambda x: bool( + ContinuousNDSupport(intervals=[Interval1D(0, np.inf)]).contains(np.array([x])) + ) + ) + family = ContinuousExponentialClassFamily( + name="ExponentialFamily", + log_partition=exponential_log_partition, + sufficient_statistics=lambda x: x, + normalization_constant=lambda _: 1, + parameter_space=support_neg, + sufficient_statistics_values=support_pos, + support=support_pos, + distr_type=UnivariateContinuous, + distr_parametrizations=["theta"], + distr_characteristics={CharacteristicName.CDF: lambda _params, x: x / (1 + x)}, + ) + + assert CharacteristicName.CDF in family.distr_characteristics + assert CharacteristicName.PDF in family.distr_characteristics + assert CharacteristicName.MEAN in family.distr_characteristics + assert CharacteristicName.VAR in family.distr_characteristics + + +def test_log_density_is_minus_infinity_outside_support( + exponential_family: ContinuousExponentialClassFamily, +) -> None: + params = ExponentialFamilyParametrization(theta=np.array([-2.0])) + log_density_func = cast( + Callable[[ExponentialFamilyParametrization, NumericArray], Number], + exponential_family.log_density, + ) + density_func = cast( + Callable[[ExponentialFamilyParametrization, NumericArray], Number], + exponential_family.density, + ) + + assert log_density_func(params, np.asarray(-0.1)) == -np.inf + assert density_func(params, np.asarray(-0.1)) == 0.0 + + +def test_distribution_rejects_theta_outside_parameter_space( + exponential_family: ContinuousExponentialClassFamily, +) -> None: + with pytest.raises(ValueError, match="theta belongs to parameter_space"): + exponential_family(theta=np.array([0.0]), parametrization_name="theta") + + +def test_transform_with_negation_moves_support_and_preserves_density( + exponential_family: ContinuousExponentialClassFamily, +) -> None: + transformed = exponential_family.transform(lambda x: -x) + params = ExponentialFamilyParametrization(theta=np.array([-1.5])) + transformed_log_density = cast( + Callable[[ExponentialFamilyParametrization, NumericArray], Number], + transformed.log_density, + ) + log_density = cast( + Callable[[ExponentialFamilyParametrization, NumericArray], Number], + exponential_family.log_density, + ) + + assert transformed.name == "TransformedExponentialFamily" + assert transformed_log_density(params, np.asarray(-2.0)) == pytest.approx( + log_density(params, np.asarray(2.0)) + ) + assert transformed_log_density(params, np.asarray(2.0)) == -np.inf + + +def test_transform_evaluates_base_measure_at_inverse_image() -> None: + family = _make_exponential_family( + normalization_constant=lambda x: cast(Number, np.asarray(x).item()) + ) + transformed = family.transform(lambda y: 2 * y) + params = ExponentialFamilyParametrization(theta=np.array([-1.5])) + transformed_log_density = cast( + Callable[[ExponentialFamilyParametrization, NumericArray], Number], + transformed.log_density, + ) + log_density = cast( + Callable[[ExponentialFamilyParametrization, NumericArray], Number], + family.log_density, + ) + + assert transformed_log_density(params, np.asarray(1.0)) == pytest.approx( + log_density(params, np.asarray(2.0)) + np.log(2.0) + ) + + +def test_posterior_hyperparameters_updates_sample_without_mutating_input( + exponential_family: ContinuousExponentialClassFamily, +) -> None: + prior = ExponentialConjugateHyperparameters( + effective_suff_stat_value=np.array([3.0]), + effective_sample_size=2.0, + ) + + posterior = exponential_family.posterior_hyperparameters(prior, sample=[0.5, 1.5]) + + assert_allclose(posterior.effective_suff_stat_value, np.array([5.0])) + assert posterior.effective_sample_size == 4.0 + assert_allclose(prior.effective_suff_stat_value, np.array([3.0])) + assert prior.effective_sample_size == 2.0 + + +def test_posterior_hyperparameters_accepts_single_observation( + exponential_family: ContinuousExponentialClassFamily, +) -> None: + prior = ExponentialConjugateHyperparameters( + effective_suff_stat_value=np.array([3.0]), + effective_sample_size=2.0, + ) + + posterior = exponential_family.posterior_hyperparameters(prior, sample=0.5) + assert_allclose(posterior.effective_suff_stat_value, np.array([3.5])) + assert posterior.effective_sample_size == 3.0 + + +def test_posterior_predictive_builds_family_with_posterior_parametrization( + exponential_family: ContinuousExponentialClassFamily, +) -> None: + predictive_family = exponential_family.posterior_predictive + + assert predictive_family.name == "PosteriorPredictiveExponentialFamily" + assert predictive_family.parametrization_names == ["posterior"] + assert predictive_family.get_parametrization("posterior") is ExponentialConjugateHyperparameters + assert CharacteristicName.PDF in predictive_family.distr_characteristics + + +@pytest.mark.parametrize( + ("xi", "nu"), + itertools.product((2.0, 3.0, 4.0), (2.0, 3.0, 4.0)), +) +def test_posterior_predictive_matches_lomax_density( + exponential_family: ContinuousExponentialClassFamily, + xi: float, + nu: float, +) -> None: + predictive = exponential_family.posterior_predictive.distribution( + parametrization_name="posterior", + effective_suff_stat_value=np.array([xi]), + effective_sample_size=nu, + ) + pdf = predictive.computation_strategy.query_method("pdf", distr=predictive) + x_values = np.array([0.25, 0.5, 1.5, 3.0, 6.0]) + + actual = np.asarray([pdf(x) for x in x_values], dtype=float).reshape(-1) + expected = np.asarray([lomax_pdf(shape=nu + 1, scale=xi, x=x) for x in x_values]) + + assert_allclose(actual, expected, rtol=1e-6) + + +def test_posterior_predictive_preserves_observation_support( + exponential_family: ContinuousExponentialClassFamily, +) -> None: + predictive = exponential_family.posterior_predictive.distribution( + parametrization_name="posterior", + effective_suff_stat_value=np.array([3.0]), + effective_sample_size=2.0, + ) + pdf = predictive.computation_strategy.query_method("pdf", distr=predictive) + + assert predictive.support is not None + assert not predictive.support.contains(np.array(-0.5)) + assert pdf(-0.5) == 0.0 + + +@pytest.mark.parametrize( + ("theta1", "theta2"), + itertools.product(range(2, 5), range(2, 5)), +) +def test_exponential_pdf(theta1, theta2, conjugate_for_exponential): + gamma_family: ContinuousExponentialClassFamily = conjugate_for_exponential + + alpha = theta2 + 1 + beta = theta1 + + exponential = gamma_family(theta=np.array([theta1, theta2]), parametrization_name="theta") + pdf = exponential.computation_strategy.query_method("pdf", distr=exponential) + + x = [i / 10 for i in range(100)] + + assert_allclose([pdf(xx) for xx in x], [gamma_pdf(alpha, beta, xx) for xx in x], rtol=1e-6) + + +@pytest.mark.parametrize( + ("theta1", "theta2"), + itertools.product(range(2, 5), range(2, 5)), +) +def test_exponential_mean(theta1, theta2, conjugate_for_exponential): + gamma_family: ContinuousExponentialClassFamily = conjugate_for_exponential + + alpha = theta2 + 1 + beta = theta1 + + exponential = gamma_family(theta=np.array([theta1, theta2]), parametrization_name="theta") + mean = exponential.computation_strategy.query_method(CharacteristicName.MEAN, distr=exponential) + assert np.isclose(mean(), alpha / beta, rtol=1e-6) + + +@pytest.mark.parametrize( + ("theta1", "theta2"), + itertools.product(range(2, 5), range(2, 5)), +) +def test_exponential_var(theta1, theta2, conjugate_for_exponential): + gamma_family: ContinuousExponentialClassFamily = conjugate_for_exponential + + alpha = theta2 + 1 + beta = theta1 + + exponential = gamma_family(theta=np.array([theta1, theta2]), parametrization_name="theta") + var = exponential.computation_strategy.query_method(CharacteristicName.VAR, distr=exponential) + assert np.isclose(var(), alpha / beta**2, rtol=1e-6)