From 1280bdf01e1c8e06909626229c1292bb2527dd6f Mon Sep 17 00:00:00 2001 From: amellgo Date: Mon, 12 May 2025 21:43:24 +0300 Subject: [PATCH 1/2] feat: Implement dynamic integrator selection, enabling the use of custom integrators --- .../g_estimation_given_mu_rqmc_based.py | 19 ++-- .../support_algorithms/integrator.py | 17 ++++ .../support_algorithms/quad_integrator.py | 27 ++++++ src/algorithms/support_algorithms/rqmc.py | 28 +++++- src/mixtures/nm_mixture.py | 90 +++++++++++-------- src/mixtures/nmv_mixture.py | 66 ++++++++------ src/mixtures/nv_mixture.py | 26 +++--- 7 files changed, 187 insertions(+), 86 deletions(-) create mode 100644 src/algorithms/support_algorithms/integrator.py create mode 100644 src/algorithms/support_algorithms/quad_integrator.py diff --git a/src/algorithms/semiparam_algorithms/nvm_semi_param_algorithms/g_estimation_given_mu_rqmc_based.py b/src/algorithms/semiparam_algorithms/nvm_semi_param_algorithms/g_estimation_given_mu_rqmc_based.py index 2f3da81..e10ba86 100644 --- a/src/algorithms/semiparam_algorithms/nvm_semi_param_algorithms/g_estimation_given_mu_rqmc_based.py +++ b/src/algorithms/semiparam_algorithms/nvm_semi_param_algorithms/g_estimation_given_mu_rqmc_based.py @@ -7,7 +7,8 @@ from scipy.integrate import quad_vec from scipy.special import gamma -from src.algorithms.support_algorithms.rqmc import RQMC +from src.algorithms.support_algorithms.integrator import Integrator +from src.algorithms.support_algorithms.rqmc import RQMCIntegrator from src.estimators.estimate_result import EstimateResult MU_DEFAULT_VALUE = 1.0 @@ -30,6 +31,7 @@ def v_sequence_default_value(n: float) -> float: GRID_SIZE_DEFAULT_VALUE: int = 200 INTEGRATION_TOLERANCE_DEFAULT_VALUE: float = 1e-2 INTEGRATION_LIMIT_DEFAULT_VALUE: int = 50 +INTEGRATOR_DEFAULT = RQMCIntegrator() class SemiParametricGEstimationGivenMuRQMCBased: @@ -69,6 +71,7 @@ def __init__(self, sample: Optional[_typing.NDArray[np.float64]] = None, **kwarg self.grid_size, self.integration_tolerance, self.integration_limit, + self.integrator ) = self._validate_kwargs(self.n, **kwargs) self.denominator: float = 2 * math.pi * self.n self.precompute_gamma_grid() @@ -78,7 +81,7 @@ def __init__(self, sample: Optional[_typing.NDArray[np.float64]] = None, **kwarg @staticmethod def _validate_kwargs( n: int, **kwargs: Unpack[ParamsAnnotation] - ) -> tuple[float, float, float, float, List[float], int, float, int]: + ) -> tuple[float, float, float, float, List[float], int, float, int, Integrator]: mu: float = kwargs.get("mu", MU_DEFAULT_VALUE) gmm: float = kwargs.get("gmm", GAMMA_DEFAULT_VALUE) u_value: float = kwargs.get("u_value", U_SEQUENCE_DEFAULT_VALUE(n)) @@ -87,7 +90,8 @@ def _validate_kwargs( grid_size: int = kwargs.get("grid_size", GRID_SIZE_DEFAULT_VALUE) integration_tolerance: float = kwargs.get("integration_tolerance", INTEGRATION_TOLERANCE_DEFAULT_VALUE) integration_limit: int = kwargs.get("integration_limit", INTEGRATION_LIMIT_DEFAULT_VALUE) - return mu, gmm, u_value, v_value, x_data, grid_size, integration_tolerance, integration_limit + integrator: Integrator = kwargs.get("integrator_type", INTEGRATOR_DEFAULT) + return mu, gmm, u_value, v_value, x_data, grid_size, integration_tolerance, integration_limit, integrator def conjugate_psi(self, u: float) -> complex: return complex((u**2) / 2, self.mu * u) @@ -162,15 +166,16 @@ def second_v_integrand(self, v: float, x: float) -> np.ndarray: x_power = self.x_powers[x][idx] return (self.second_u_integrals[idx] * x_power) / gamma_val - def compute_integrals_for_x(self, x: float) -> float: + def compute_integrals_for_x(self, x: float, integrator: Integrator = None) -> float: """Compute integrals using RQMC for v-integration.""" - first_integral = RQMC(lambda t: np.sum(self.first_v_integrand(t * self.v_value, x)) * self.v_value).rqmc()[0] + integrator = integrator or self.integrator + first_integral = integrator.compute_integral(func=lambda t: np.sum(self.first_v_integrand(t * self.v_value, x)) * self.v_value).value - second_integral = RQMC(lambda t: np.sum(self.second_v_integrand(-t * self.v_value, x)) * self.v_value).rqmc()[0] + second_integral = integrator.compute_integral(func=lambda t: np.sum(self.second_v_integrand(-t * self.v_value, x)) * self.v_value).value total = (first_integral + second_integral) / self.denominator return max(0.0, total.real) def algorithm(self, sample: np._typing.NDArray) -> EstimateResult: - y_data = [self.compute_integrals_for_x(x) for x in self.x_data] + y_data = [self.compute_integrals_for_x(x, self.integrator) for x in self.x_data] return EstimateResult(list_value=y_data, success=True) diff --git a/src/algorithms/support_algorithms/integrator.py b/src/algorithms/support_algorithms/integrator.py new file mode 100644 index 0000000..6731aed --- /dev/null +++ b/src/algorithms/support_algorithms/integrator.py @@ -0,0 +1,17 @@ +from dataclasses import dataclass +from typing import Any, Protocol, Callable, Optional + + +@dataclass +class IntegrationResult: + value: float + error: float + message: Optional[dict[str, Any]] | None = None + + +class Integrator(Protocol): + + """Base class for integral calculation""" + + def compute_integral(self, func: Callable, params: dict) -> IntegrationResult: + ... diff --git a/src/algorithms/support_algorithms/quad_integrator.py b/src/algorithms/support_algorithms/quad_integrator.py new file mode 100644 index 0000000..937cb46 --- /dev/null +++ b/src/algorithms/support_algorithms/quad_integrator.py @@ -0,0 +1,27 @@ +from typing import Callable + +from scipy.integrate import quad +from src.algorithms.support_algorithms.integrator import IntegrationResult + +class QuadIntegrator: + + def compute_integral(self, func: Callable, params: dict) -> IntegrationResult: + + """ + Compute integral via quad integrator + + Args: + func: integrated function + params: Parameters of integration algorithm + + Returns: moment approximation and error tolerance + """ + + full_output_requested = params.pop('full_output', False) + quad_res = quad(func, **params) + if full_output_requested: + value, error, message = quad_res + else: + value, error = quad_res + message = None + return IntegrationResult(value, error, message) diff --git a/src/algorithms/support_algorithms/rqmc.py b/src/algorithms/support_algorithms/rqmc.py index 4b10c9d..520ba62 100644 --- a/src/algorithms/support_algorithms/rqmc.py +++ b/src/algorithms/support_algorithms/rqmc.py @@ -5,6 +5,8 @@ import scipy from numba import njit +from src.algorithms.support_algorithms.integrator import IntegrationResult + BITS = 30 """Number of bits in XOR. Should be less than 64""" NUMBA_FAST_MATH = True @@ -126,7 +128,8 @@ def _update( Returns:Updated mean of all rows - """ + + """ values = [] sum_of_new: float = 0.0 for i in range(self.count): @@ -212,9 +215,9 @@ def _xor_float(a: float, b: float) -> float: Returns: XOR float value """ - a = int(a * (2**BITS)) - b = int(b * (2**BITS)) - return np.bitwise_xor(a, b) / 2**BITS + a = int(a * (2 ** BITS)) + b = int(b * (2 ** BITS)) + return np.bitwise_xor(a, b) / 2 ** BITS def __call__(self) -> tuple[float, float]: """Interface for users @@ -223,3 +226,20 @@ def __call__(self) -> tuple[float, float]: """ return self.rqmc() + + +class RQMCIntegrator: + + def compute_integral(self, func: Callable, params: dict) -> IntegrationResult: + """ + Compute integral via RQMC integrator + + Args: + func: integrated function + params: Parameters of integration algorithm + + Returns: moment approximation and error tolerance + """ + + rqmc = RQMC(func, **params)() + return IntegrationResult(rqmc[0], rqmc[1]) diff --git a/src/mixtures/nm_mixture.py b/src/mixtures/nm_mixture.py index 33c7580..0e4d811 100644 --- a/src/mixtures/nm_mixture.py +++ b/src/mixtures/nm_mixture.py @@ -2,13 +2,14 @@ from typing import Any import numpy as np -from scipy.integrate import quad from scipy.special import binom from scipy.stats import norm, rv_continuous from scipy.stats.distributions import rv_frozen from src.algorithms.support_algorithms.log_rqmc import LogRQMC -from src.algorithms.support_algorithms.rqmc import RQMC +from src.algorithms.support_algorithms.integrator import Integrator +from src.algorithms.support_algorithms.rqmc import RQMCIntegrator +from src.algorithms.support_algorithms.quad_integrator import QuadIntegrator from src.mixtures.abstract_mixture import AbstractMixtures @@ -59,156 +60,175 @@ def _params_validation(self, data_collector: Any, params: dict[str, float | rv_c raise ValueError("Gamma cant be zero") return data_class - def _classical_moment(self, n: int, params: dict) -> tuple[float, float]: + def _classical_moment(self, n: int, params: dict, integrator: Integrator = None) -> tuple[float, float]: """ Compute n-th moment of classical NMM Args: n (): Moment ordinal params (): Parameters of integration algorithm + integrator (): type of integrator to computing Returns: moment approximation and error tolerance """ mixture_moment = 0 error_tolerance = 0 + integrator = integrator or QuadIntegrator() for k in range(0, n + 1): for l in range(0, k + 1): - coefficient = binom(n, n - k) * binom(k, k - l) * (self.params.beta ** (k - l)) * (self.params.gamma**l) - mixing_moment = quad(lambda u: self.params.distribution.ppf(u) ** (k - l), 0, 1, **params) - error_tolerance += (self.params.beta ** (k - l)) * mixing_moment[1] - mixture_moment += coefficient * (self.params.alpha ** (n - k)) * mixing_moment[0] * norm.moment(l) - return mixture_moment, error_tolerance + coefficient = binom(n, n - k) * binom(k, k - l) * (self.params.beta ** (k - l)) * ( + self.params.gamma ** l) + mixing_moment = integrator.compute_integral(lambda u: self.params.distribution.ppf(u) ** (k - l), + {"a": 0, "b": 1, **params}) + error_tolerance += (self.params.beta ** (k - l)) * mixing_moment.error + mixture_moment += coefficient * (self.params.alpha ** (n - k)) * mixing_moment.value * norm.moment(l) + return mixture_moment, error_tolerance + - def _canonical_moment(self, n: int, params: dict) -> tuple[float, float]: + def _canonical_moment(self, n: int, params: dict, integrator: Integrator = None) -> tuple[float, float]: """ Compute n-th moment of canonical NMM Args: n (): Moment ordinal params (): Parameters of integration algorithm + integrator (): type of integrator to computing Returns: moment approximation and error tolerance """ mixture_moment = 0 error_tolerance = 0 + integrator = integrator or QuadIntegrator() for k in range(0, n + 1): - coefficient = binom(n, n - k) * (self.params.sigma**k) - mixing_moment = quad(lambda u: self.params.distribution.ppf(u) ** (n - k), 0, 1, **params) - error_tolerance += mixing_moment[1] - mixture_moment += coefficient * mixing_moment[0] * norm.moment(k) + coefficient = binom(n, n - k) * (self.params.sigma ** k) + mixing_moment = integrator.compute_integral(lambda u: self.params.distribution.ppf(u) ** (n - k), + {"a": 0, "b": 1, **params}) + error_tolerance += mixing_moment.error + mixture_moment += coefficient * mixing_moment.value * norm.moment(k) return mixture_moment, error_tolerance - def compute_moment(self, n: int, params: dict) -> tuple[float, float]: + def compute_moment(self, n: int, params: dict, integrator: Integrator = None) -> tuple[float, float]: """ Compute n-th moment of NMM Args: n (): Moment ordinal params (): Parameters of integration algorithm + integrator (): type of integrator to computing Returns: moment approximation and error tolerance """ if isinstance(self.params, _NMMClassicDataCollector): - return self._classical_moment(n, params) - return self._canonical_moment(n, params) + return self._classical_moment(n, params, integrator) + return self._canonical_moment(n, params, integrator) - def _canonical_compute_cdf(self, x: float, params: dict) -> tuple[float, float]: + def _canonical_compute_cdf(self, x: float, params: dict, integrator: Integrator = None) -> tuple[float, float]: """ Equation for canonical cdf Args: x (): point params (): parameters of RQMC algorithm + integrator (): type of integrator to computing Returns: computed cdf and error tolerance """ - rqmc = RQMC(lambda u: norm.cdf((x - self.params.distribution.ppf(u)) / np.abs(self.params.sigma)), **params) - return rqmc() + integrator = integrator or RQMCIntegrator() + rqmc = integrator.compute_integral(func=lambda u: norm.cdf((x - self.params.distribution.ppf(u)) / np.abs(self.params.sigma)), **params) + return rqmc.value, rqmc.error - def _classical_compute_cdf(self, x: float, params: dict) -> tuple[float, float]: + def _classical_compute_cdf(self, x: float, params: dict, integrator: Integrator = None) -> tuple[float, float]: """ Equation for classic cdf Args: x (): point params (): parameters of RQMC algorithm + integrator (): type of integrator to computing Returns: computed cdf and error tolerance """ - rqmc = RQMC( + integrator = integrator or RQMCIntegrator() + rqmc = integrator.compute_integral(func= lambda u: norm.cdf( (x - self.params.alpha - self.params.beta * self.params.distribution.ppf(u)) / np.abs(self.params.gamma) ), **params ) - return rqmc() + return rqmc.value, rqmc.error - def compute_cdf(self, x: float, params: dict) -> tuple[float, float]: + def compute_cdf(self, x: float, params: dict, integrator: Integrator = None) -> tuple[float, float]: """ Choose equation for cdf estimation depends on Mixture form Args: x (): point params (): parameters of RQMC algorithm + integrator (): type of integrator to computing Returns: Computed pdf and error tolerance """ if isinstance(self.params, _NMMCanonicalDataCollector): - return self._canonical_compute_cdf(x, params) - return self._classical_compute_cdf(x, params) + return self._canonical_compute_cdf(x, params, integrator) + return self._classical_compute_cdf(x, params, integrator) - def _canonical_compute_pdf(self, x: float, params: dict) -> tuple[float, float]: + def _canonical_compute_pdf(self, x: float, params: dict, integrator: Integrator = None) -> tuple[float, float]: """ Equation for canonical pdf Args: x (): point params (): parameters of RQMC algorithm + integrator (): type of integrator to computing Returns: computed pdf and error tolerance """ - rqmc = RQMC( + integrator = integrator or RQMCIntegrator() + rqmc = integrator.compute_integral(func= lambda u: (1 / np.abs(self.params.sigma)) * norm.pdf((x - self.params.distribution.ppf(u)) / np.abs(self.params.sigma)), **params ) - return rqmc() + return rqmc.value, rqmc.error - def _classical_compute_pdf(self, x: float, params: dict) -> tuple[float, float]: + def _classical_compute_pdf(self, x: float, params: dict, integrator: Integrator = None) -> tuple[float, float]: """ Equation for classic pdf Args: x (): point params (): parameters of RQMC algorithm + integrator (): type of integrator to computing Returns: computed pdf and error tolerance """ - rqmc = RQMC( + integrator = integrator or RQMCIntegrator() + rqmc = integrator.compute_integral(func= lambda u: (1 / np.abs(self.params.gamma)) * norm.pdf( (x - self.params.alpha - self.params.beta * self.params.distribution.ppf(u)) / np.abs(self.params.gamma) ), **params ) - return rqmc() + return rqmc.value, rqmc.error - def compute_pdf(self, x: float, params: dict) -> tuple[float, float]: + def compute_pdf(self, x: float, params: dict, integrator: Integrator = None) -> tuple[float, float]: """ Choose equation for pdf estimation depends on Mixture form Args: x (): point params (): parameters of RQMC algorithm + integrator (): type of integrator to computing Returns: Computed pdf and error tolerance """ if isinstance(self.params, _NMMCanonicalDataCollector): - return self._canonical_compute_pdf(x, params) - return self._classical_compute_pdf(x, params) + return self._canonical_compute_pdf(x, params, integrator) + return self._classical_compute_pdf(x, params, integrator) def _classical_compute_log_pdf(self, x: float, params: dict) -> tuple[float, float]: """ @@ -258,4 +278,4 @@ def compute_logpdf(self, x: float, params: dict) -> tuple[float, float]: """ if isinstance(self.params, _NMMCanonicalDataCollector): return self._canonical_compute_log_pdf(x, params) - return self._classical_compute_log_pdf(x, params) + return self._classical_compute_log_pdf(x, params) \ No newline at end of file diff --git a/src/mixtures/nmv_mixture.py b/src/mixtures/nmv_mixture.py index c090ff3..9f2a7e0 100644 --- a/src/mixtures/nmv_mixture.py +++ b/src/mixtures/nmv_mixture.py @@ -8,7 +8,8 @@ from scipy.stats.distributions import rv_frozen from src.algorithms.support_algorithms.log_rqmc import LogRQMC -from src.algorithms.support_algorithms.rqmc import RQMC +from src.algorithms.support_algorithms.integrator import Integrator +from src.algorithms.support_algorithms.rqmc import RQMCIntegrator from src.mixtures.abstract_mixture import AbstractMixtures @@ -40,13 +41,14 @@ class NormalMeanVarianceMixtures(AbstractMixtures): def __init__(self, mixture_form: str, **kwargs: Any) -> None: super().__init__(mixture_form, **kwargs) - def _classical_moment(self, n: int, params: dict) -> tuple[float, float]: + def _classical_moment(self, n: int, params: dict, integrator: Integrator = None) -> tuple[float, float]: """ Compute n-th moment of classical NMM Args: n (): Moment ordinal params (): Parameters of integration algorithm + integrator (): type of integrator to computing Returns: moment approximation and error tolerance @@ -67,16 +69,18 @@ def integral_func(u: float) -> float: ) return result - rqmc = RQMC(lambda u: integral_func(u), **params) - return rqmc() + integrator = integrator or RQMCIntegrator() + rqmc = integrator.compute_integral(func=lambda u: integral_func(u), **params) + return rqmc.value, rqmc.error - def _canonical_moment(self, n: int, params: dict) -> tuple[float, float]: + def _canonical_moment(self, n: int, params: dict, integrator: Integrator = None) -> tuple[float, float]: """ Compute n-th moment of classical NMM Args: n (): Moment ordinal params (): Parameters of integration algorithm + integrator (): type of integrator to computing Returns: moment approximation and error tolerance @@ -96,15 +100,16 @@ def integral_func(u: float) -> float: ) return result - rqmc = RQMC(lambda u: integral_func(u), **params) - return rqmc() + integrator = integrator or RQMCIntegrator() + rqmc = integrator.compute_integral(func=lambda u: integral_func(u), **params) + return rqmc.value, rqmc.error - def compute_moment(self, n: int, params: dict) -> tuple[float, float]: + def compute_moment(self, n: int, params: dict, integrator: Integrator = None) -> tuple[float, float]: if isinstance(self.params, _NMVMClassicDataCollector): - return self._classical_moment(n, params) - return self._canonical_moment(n, params) + return self._classical_moment(n, params, integrator) + return self._canonical_moment(n, params, integrator) - def _classical_cdf(self, x: float, params: dict) -> tuple[float, float]: + def _classical_cdf(self, x: float, params: dict, integrator: Integrator = None) -> tuple[float, float]: def _inner_func(u: float) -> float: ppf = lru_cache()(self.params.distribution.ppf)(u) point = (x - self.params.alpha) / (np.sqrt(ppf) * self.params.gamma) - ( @@ -112,24 +117,26 @@ def _inner_func(u: float) -> float: ) return norm.cdf(point) - rqmc = RQMC(lambda u: _inner_func(u), **params) - return rqmc() + integrator = integrator or RQMCIntegrator() + rqmc = integrator.compute_integral(func=lambda u: _inner_func(u), **params) + return rqmc.value, rqmc.error - def _canonical_cdf(self, x: float, params: dict) -> tuple[float, float]: + def _canonical_cdf(self, x: float, params: dict, integrator: Integrator = None) -> tuple[float, float]: def _inner_func(u: float) -> float: ppf = self.params.distribution.ppf(u) point = (x - self.params.alpha) / (np.sqrt(ppf)) - (self.params.mu * np.sqrt(ppf)) return norm.cdf(point) - rqmc = RQMC(lambda u: _inner_func(u), **params) - return rqmc() + integrator = integrator or RQMCIntegrator() + rqmc = integrator.compute_integral(func=lambda u: _inner_func(u), **params) + return rqmc.value, rqmc.error - def compute_cdf(self, x: float, params: dict) -> tuple[float, float]: + def compute_cdf(self, x: float, params: dict, integrator: Integrator = None) -> tuple[float, float]: if isinstance(self.params, _NMVMClassicDataCollector): - return self._classical_cdf(x, params) - return self._canonical_cdf(x, params) + return self._classical_cdf(x, params, integrator) + return self._canonical_cdf(x, params, integrator) - def _classical_pdf(self, x: float, params: dict) -> tuple[float, float]: + def _classical_pdf(self, x: float, params: dict, integrator: Integrator = None) -> tuple[float, float]: def _inner_func(u: float) -> float: ppf = self.params.distribution.ppf(u) return ( @@ -140,10 +147,11 @@ def _inner_func(u: float) -> float: ) ) - rqmc = RQMC(lambda u: _inner_func(u), **params)() - return np.exp(self.params.beta * (x - self.params.alpha) / self.params.gamma**2) * rqmc[0], rqmc[1] + integrator = integrator or RQMCIntegrator() + rqmc = integrator.compute_integral(func=lambda u: _inner_func(u), **params) + return np.exp(self.params.beta * (x - self.params.alpha) / self.params.gamma**2) * rqmc.value, rqmc.error - def _canonical_pdf(self, x: float, params: dict) -> tuple[float, float]: + def _canonical_pdf(self, x: float, params: dict, integrator: Integrator = None) -> tuple[float, float]: def _inner_func(u: float) -> float: ppf = self.params.distribution.ppf(u) return ( @@ -152,9 +160,9 @@ def _inner_func(u: float) -> float: * np.exp(-((x - self.params.alpha) ** 2 + self.params.mu**2 * ppf**2) / (2 * ppf)) ) - rqmc = RQMC(lambda u: _inner_func(u), **params) - res = rqmc() - return np.exp(self.params.mu * (x - self.params.alpha)) * res[0], res[1] + integrator = integrator or RQMCIntegrator() + rqmc = integrator.compute_integral(func=lambda u: _inner_func(u), **params) + return np.exp(self.params.mu * (x - self.params.alpha)) * rqmc.value, rqmc.error def _classical_log_pdf(self, x: float, params: dict) -> tuple[float, float]: def _inner_func(u: float) -> float: @@ -178,10 +186,10 @@ def _inner_func(u: float) -> float: rqmc = LogRQMC(lambda u: _inner_func(u), **params) return rqmc() - def compute_pdf(self, x: float, params: dict) -> tuple[float, float]: + def compute_pdf(self, x: float, params: dict, integrator: Integrator = None) -> tuple[float, float]: if isinstance(self.params, _NMVMClassicDataCollector): - return self._classical_pdf(x, params) - return self._canonical_pdf(x, params) + return self._classical_pdf(x, params, integrator) + return self._canonical_pdf(x, params, integrator) def compute_logpdf(self, x: float, params: dict) -> tuple[Any, float]: if isinstance(self.params, _NMVMClassicDataCollector): diff --git a/src/mixtures/nv_mixture.py b/src/mixtures/nv_mixture.py index 7905e6e..1be3343 100644 --- a/src/mixtures/nv_mixture.py +++ b/src/mixtures/nv_mixture.py @@ -8,7 +8,8 @@ from scipy.stats.distributions import rv_frozen from src.algorithms.support_algorithms.log_rqmc import LogRQMC -from src.algorithms.support_algorithms.rqmc import RQMC +from src.algorithms.support_algorithms.integrator import Integrator +from src.algorithms.support_algorithms.rqmc import RQMCIntegrator from src.mixtures.abstract_mixture import AbstractMixtures @@ -39,12 +40,13 @@ class NormalVarianceMixtures(AbstractMixtures): def __init__(self, mixture_form: str, **kwargs: Any) -> None: super().__init__(mixture_form, **kwargs) - def compute_moment(self, n: int, params: dict) -> tuple[float, float]: + def compute_moment(self, n: int, params: dict, integrator: Integrator = None) -> tuple[float, float]: """ Compute n-th moment of NVM Args: n (): Moment ordinal params (): Parameters of integration algorithm + integrator (): type of integrator to computing Returns: moment approximation and error tolerance """ gamma = self.params.gamma if isinstance(self.params, _NVMClassicDataCollector) else 1 @@ -60,16 +62,17 @@ def integrate_func(u: float) -> float: for k in range(0, n + 1) ] ) + integrator = integrator or RQMCIntegrator() + result = integrator.compute_integral(func=integrate_func, **params) + return result.value, result.error - result = RQMC(integrate_func, **params)() - return result - - def compute_cdf(self, x: float, params: dict) -> tuple[float, float]: + def compute_cdf(self, x: float, params: dict, integrator: Integrator = None) -> tuple[float, float]: parametric_norm = norm(0, self.params.gamma if isinstance(self.params, _NVMClassicDataCollector) else 1) - rqmc = RQMC( + integrator = integrator or RQMCIntegrator() + rqmc = integrator.compute_integral(func= lambda u: parametric_norm.cdf((x - self.params.alpha) / np.sqrt(self.params.distribution.ppf(u))), **params ) - return rqmc() + return rqmc.value, rqmc.error @lru_cache() def _integrand_func(self, u: float, d: float, gamma: float) -> float: @@ -80,11 +83,12 @@ def _log_integrand_func(self, u: float, d: float, gamma: float | int | np.int64) ppf = self.params.distribution.ppf(u) return -(ppf * np.log(np.pi * 2 * ppf * gamma**2) + d) / (2 * ppf) - def compute_pdf(self, x: float, params: dict) -> tuple[float, float]: + def compute_pdf(self, x: float, params: dict, integrator: Integrator = None) -> tuple[float, float]: gamma = self.params.gamma if isinstance(self.params, _NVMClassicDataCollector) else 1 d = (x - self.params.alpha) ** 2 / gamma**2 - rqmc = RQMC(lambda u: self._integrand_func(u, d, gamma), **params) - return rqmc() + integrator = integrator or RQMCIntegrator() + rqmc = integrator.compute_integral(func=lambda u: self._integrand_func(u, d, gamma), **params) + return rqmc.value, rqmc.error def compute_logpdf(self, x: float, params: dict) -> tuple[float, float]: gamma = self.params.gamma if isinstance(self.params, _NVMClassicDataCollector) else 1 From 51a7303c04ef071f415132105d8e3176d1329d4d Mon Sep 17 00:00:00 2001 From: amellgo Date: Sun, 1 Jun 2025 23:10:54 +0300 Subject: [PATCH 2/2] feat: integration parameters are now part of the constructor --- .../g_estimation_given_mu_rqmc_based.py | 18 ++-- .../support_algorithms/integrator.py | 5 +- .../support_algorithms/quad_integrator.py | 49 +++++++-- src/algorithms/support_algorithms/rqmc.py | 39 ++++++- src/mixtures/abstract_mixture.py | 7 +- src/mixtures/nm_mixture.py | 101 +++++++----------- src/mixtures/nmv_mixture.py | 74 ++++++------- src/mixtures/nv_mixture.py | 24 ++--- 8 files changed, 174 insertions(+), 143 deletions(-) diff --git a/src/algorithms/semiparam_algorithms/nvm_semi_param_algorithms/g_estimation_given_mu_rqmc_based.py b/src/algorithms/semiparam_algorithms/nvm_semi_param_algorithms/g_estimation_given_mu_rqmc_based.py index e10ba86..fc089cf 100644 --- a/src/algorithms/semiparam_algorithms/nvm_semi_param_algorithms/g_estimation_given_mu_rqmc_based.py +++ b/src/algorithms/semiparam_algorithms/nvm_semi_param_algorithms/g_estimation_given_mu_rqmc_based.py @@ -31,7 +31,6 @@ def v_sequence_default_value(n: float) -> float: GRID_SIZE_DEFAULT_VALUE: int = 200 INTEGRATION_TOLERANCE_DEFAULT_VALUE: float = 1e-2 INTEGRATION_LIMIT_DEFAULT_VALUE: int = 50 -INTEGRATOR_DEFAULT = RQMCIntegrator() class SemiParametricGEstimationGivenMuRQMCBased: @@ -70,8 +69,7 @@ def __init__(self, sample: Optional[_typing.NDArray[np.float64]] = None, **kwarg self.x_data, self.grid_size, self.integration_tolerance, - self.integration_limit, - self.integrator + self.integration_limit ) = self._validate_kwargs(self.n, **kwargs) self.denominator: float = 2 * math.pi * self.n self.precompute_gamma_grid() @@ -81,7 +79,7 @@ def __init__(self, sample: Optional[_typing.NDArray[np.float64]] = None, **kwarg @staticmethod def _validate_kwargs( n: int, **kwargs: Unpack[ParamsAnnotation] - ) -> tuple[float, float, float, float, List[float], int, float, int, Integrator]: + ) -> tuple[float, float, float, float, List[float], int, float, int]: mu: float = kwargs.get("mu", MU_DEFAULT_VALUE) gmm: float = kwargs.get("gmm", GAMMA_DEFAULT_VALUE) u_value: float = kwargs.get("u_value", U_SEQUENCE_DEFAULT_VALUE(n)) @@ -90,8 +88,7 @@ def _validate_kwargs( grid_size: int = kwargs.get("grid_size", GRID_SIZE_DEFAULT_VALUE) integration_tolerance: float = kwargs.get("integration_tolerance", INTEGRATION_TOLERANCE_DEFAULT_VALUE) integration_limit: int = kwargs.get("integration_limit", INTEGRATION_LIMIT_DEFAULT_VALUE) - integrator: Integrator = kwargs.get("integrator_type", INTEGRATOR_DEFAULT) - return mu, gmm, u_value, v_value, x_data, grid_size, integration_tolerance, integration_limit, integrator + return mu, gmm, u_value, v_value, x_data, grid_size, integration_tolerance, integration_limit def conjugate_psi(self, u: float) -> complex: return complex((u**2) / 2, self.mu * u) @@ -166,16 +163,15 @@ def second_v_integrand(self, v: float, x: float) -> np.ndarray: x_power = self.x_powers[x][idx] return (self.second_u_integrals[idx] * x_power) / gamma_val - def compute_integrals_for_x(self, x: float, integrator: Integrator = None) -> float: + def compute_integrals_for_x(self, x: float, integrator: Integrator = RQMCIntegrator()) -> float: """Compute integrals using RQMC for v-integration.""" - integrator = integrator or self.integrator - first_integral = integrator.compute_integral(func=lambda t: np.sum(self.first_v_integrand(t * self.v_value, x)) * self.v_value).value + first_integral = integrator.compute(func=lambda t: np.sum(self.first_v_integrand(t * self.v_value, x)) * self.v_value).value - second_integral = integrator.compute_integral(func=lambda t: np.sum(self.second_v_integrand(-t * self.v_value, x)) * self.v_value).value + second_integral = integrator.compute(func=lambda t: np.sum(self.second_v_integrand(-t * self.v_value, x)) * self.v_value).value total = (first_integral + second_integral) / self.denominator return max(0.0, total.real) def algorithm(self, sample: np._typing.NDArray) -> EstimateResult: - y_data = [self.compute_integrals_for_x(x, self.integrator) for x in self.x_data] + y_data = [self.compute_integrals_for_x(x, RQMCIntegrator()) for x in self.x_data] return EstimateResult(list_value=y_data, success=True) diff --git a/src/algorithms/support_algorithms/integrator.py b/src/algorithms/support_algorithms/integrator.py index 6731aed..8ee0fb7 100644 --- a/src/algorithms/support_algorithms/integrator.py +++ b/src/algorithms/support_algorithms/integrator.py @@ -13,5 +13,8 @@ class Integrator(Protocol): """Base class for integral calculation""" - def compute_integral(self, func: Callable, params: dict) -> IntegrationResult: + def __init__(self) -> None: + ... + + def compute(self, func: Callable) -> IntegrationResult: ... diff --git a/src/algorithms/support_algorithms/quad_integrator.py b/src/algorithms/support_algorithms/quad_integrator.py index 937cb46..616eca2 100644 --- a/src/algorithms/support_algorithms/quad_integrator.py +++ b/src/algorithms/support_algorithms/quad_integrator.py @@ -1,27 +1,60 @@ -from typing import Callable +from typing import Callable, Any, Sequence from scipy.integrate import quad from src.algorithms.support_algorithms.integrator import IntegrationResult class QuadIntegrator: - def compute_integral(self, func: Callable, params: dict) -> IntegrationResult: + def __init__( + self, + a: float = 0, + b: float = 1, + args: tuple[Any, ...] = (), + full_output: int = 0, + epsabs: float | int = 1.49e-08, + epsrel: float | int = 1.49e-08, + limit: float | int = 50, + points: Sequence[float | int] | None = None, + weight: float | int | None = None, + wvar: Any = None, + wopts: Any = None, + maxp1: float | int = 50, + limlst: int = 50, + complex_func: bool = False, + ): + self.params = { + 'a': a, + 'b': b, + 'args': args, + 'full_output': full_output, + 'epsabs': epsabs, + 'epsrel': epsrel, + 'limit': limit, + 'points': points, + 'weight': weight, + 'wvar': wvar, + 'wopts': wopts, + 'maxp1': maxp1, + 'limlst': limlst, + 'complex_func': complex_func + } + + def compute(self, func: Callable) -> IntegrationResult: """ Compute integral via quad integrator Args: func: integrated function - params: Parameters of integration algorithm Returns: moment approximation and error tolerance """ - full_output_requested = params.pop('full_output', False) - quad_res = quad(func, **params) - if full_output_requested: - value, error, message = quad_res + verbose = self.params.pop('full_output', False) + result = quad(func, **self.params) + if verbose: + value, error, message = result else: - value, error = quad_res + value, error = result message = None return IntegrationResult(value, error, message) diff --git a/src/algorithms/support_algorithms/rqmc.py b/src/algorithms/support_algorithms/rqmc.py index 520ba62..543d2cb 100644 --- a/src/algorithms/support_algorithms/rqmc.py +++ b/src/algorithms/support_algorithms/rqmc.py @@ -229,17 +229,48 @@ def __call__(self) -> tuple[float, float]: class RQMCIntegrator: + """ + Randomize Quasi Monte Carlo Method + + Args: + error_tolerance: pre-specified error tolerance + count: number of rows of random values matrix + base_n: number of columns of random values matrix + i_max: allowed number of cycles + a: parameter for quantile of normal distribution + + """ + + def __init__( + self, + error_tolerance: float = 1e-6, + count: int = 25, + base_n: int = 2 ** 6, + i_max: int = 100, + a: float = 0.00047, + ): + self.error_tolerance = error_tolerance + self.count = count + self.base_n = base_n + self.i_max = i_max + self.a = a - def compute_integral(self, func: Callable, params: dict) -> IntegrationResult: + def compute(self, func: Callable) -> IntegrationResult: """ Compute integral via RQMC integrator Args: func: integrated function - params: Parameters of integration algorithm Returns: moment approximation and error tolerance """ + result = RQMC( + func, + error_tolerance=self.error_tolerance, + count=self.count, + base_n=self.base_n, + i_max=self.i_max, + a=self.a, + )() + return IntegrationResult(result[0], result[1]) - rqmc = RQMC(func, **params)() - return IntegrationResult(rqmc[0], rqmc[1]) diff --git a/src/mixtures/abstract_mixture.py b/src/mixtures/abstract_mixture.py index 13c8e42..3e95de4 100644 --- a/src/mixtures/abstract_mixture.py +++ b/src/mixtures/abstract_mixture.py @@ -5,6 +5,7 @@ from scipy.stats import rv_continuous from scipy.stats.distributions import rv_frozen +from src.algorithms.support_algorithms.integrator import Integrator class AbstractMixtures(metaclass=ABCMeta): """Base class for Mixtures""" @@ -28,13 +29,13 @@ def __init__(self, mixture_form: str, **kwargs: Any) -> None: raise AssertionError(f"Unknown mixture form: {mixture_form}") @abstractmethod - def compute_moment(self, n: int, params: dict) -> tuple[float, float]: ... + def compute_moment(self, n: int, integrator: Integrator) -> tuple[float, float]: ... @abstractmethod - def compute_cdf(self, x: float, params: dict) -> tuple[float, float]: ... + def compute_cdf(self, x: float, integrator: Integrator) -> tuple[float, float]: ... @abstractmethod - def compute_pdf(self, x: float, params: dict) -> tuple[float, float]: ... + def compute_pdf(self, x: float, integrator: Integrator) -> tuple[float, float]: ... @abstractmethod def compute_logpdf(self, x: float, params: dict) -> tuple[float, float]: ... diff --git a/src/mixtures/nm_mixture.py b/src/mixtures/nm_mixture.py index 0e4d811..617c498 100644 --- a/src/mixtures/nm_mixture.py +++ b/src/mixtures/nm_mixture.py @@ -60,175 +60,154 @@ def _params_validation(self, data_collector: Any, params: dict[str, float | rv_c raise ValueError("Gamma cant be zero") return data_class - def _classical_moment(self, n: int, params: dict, integrator: Integrator = None) -> tuple[float, float]: + def _classical_moment(self, n: int, integrator: Integrator = QuadIntegrator()) -> tuple[float, float]: """ Compute n-th moment of classical NMM Args: n (): Moment ordinal - params (): Parameters of integration algorithm - integrator (): type of integrator to computing + integrator (): class of integrator with params to computing Returns: moment approximation and error tolerance """ - mixture_moment = 0 - error_tolerance = 0 - integrator = integrator or QuadIntegrator() + mixture_moment = 0.0 + error_tolerance = 0.0 for k in range(0, n + 1): for l in range(0, k + 1): coefficient = binom(n, n - k) * binom(k, k - l) * (self.params.beta ** (k - l)) * ( self.params.gamma ** l) - mixing_moment = integrator.compute_integral(lambda u: self.params.distribution.ppf(u) ** (k - l), - {"a": 0, "b": 1, **params}) + mixing_moment = integrator.compute(lambda u: self.params.distribution.ppf(u) ** (k - l)) error_tolerance += (self.params.beta ** (k - l)) * mixing_moment.error mixture_moment += coefficient * (self.params.alpha ** (n - k)) * mixing_moment.value * norm.moment(l) - return mixture_moment, error_tolerance - + return mixture_moment, error_tolerance - def _canonical_moment(self, n: int, params: dict, integrator: Integrator = None) -> tuple[float, float]: + def _canonical_moment(self, n: int, integrator: Integrator = QuadIntegrator()) -> tuple[float, float]: """ Compute n-th moment of canonical NMM Args: n (): Moment ordinal - params (): Parameters of integration algorithm - integrator (): type of integrator to computing + integrator (): class of integrator with params to computing Returns: moment approximation and error tolerance """ - mixture_moment = 0 - error_tolerance = 0 - integrator = integrator or QuadIntegrator() + mixture_moment = 0.0 + error_tolerance = 0.0 for k in range(0, n + 1): coefficient = binom(n, n - k) * (self.params.sigma ** k) - mixing_moment = integrator.compute_integral(lambda u: self.params.distribution.ppf(u) ** (n - k), - {"a": 0, "b": 1, **params}) + mixing_moment = integrator.compute(lambda u: self.params.distribution.ppf(u) ** (n - k)) error_tolerance += mixing_moment.error mixture_moment += coefficient * mixing_moment.value * norm.moment(k) return mixture_moment, error_tolerance - def compute_moment(self, n: int, params: dict, integrator: Integrator = None) -> tuple[float, float]: + def compute_moment(self, n: int, integrator: Integrator = QuadIntegrator()) -> tuple[float, float]: """ Compute n-th moment of NMM Args: n (): Moment ordinal - params (): Parameters of integration algorithm - integrator (): type of integrator to computing + integrator (): class of integrator with params to computing Returns: moment approximation and error tolerance """ if isinstance(self.params, _NMMClassicDataCollector): - return self._classical_moment(n, params, integrator) - return self._canonical_moment(n, params, integrator) + return self._classical_moment(n, integrator) + return self._canonical_moment(n, integrator) - def _canonical_compute_cdf(self, x: float, params: dict, integrator: Integrator = None) -> tuple[float, float]: + def _canonical_compute_cdf(self, x: float, integrator: Integrator = RQMCIntegrator()) -> tuple[float, float]: """ Equation for canonical cdf Args: x (): point - params (): parameters of RQMC algorithm - integrator (): type of integrator to computing + integrator (): class of integrator with params to computing Returns: computed cdf and error tolerance """ - integrator = integrator or RQMCIntegrator() - rqmc = integrator.compute_integral(func=lambda u: norm.cdf((x - self.params.distribution.ppf(u)) / np.abs(self.params.sigma)), **params) - return rqmc.value, rqmc.error + result = integrator.compute(func=lambda u: norm.cdf((x - self.params.distribution.ppf(u)) / np.abs(self.params.sigma))) + return result.value, result.error - def _classical_compute_cdf(self, x: float, params: dict, integrator: Integrator = None) -> tuple[float, float]: + def _classical_compute_cdf(self, x: float, integrator: Integrator = RQMCIntegrator()) -> tuple[float, float]: """ Equation for classic cdf Args: x (): point - params (): parameters of RQMC algorithm - integrator (): type of integrator to computing + integrator (): class of integrator with params to computing Returns: computed cdf and error tolerance """ - integrator = integrator or RQMCIntegrator() - rqmc = integrator.compute_integral(func= + result = integrator.compute(func= lambda u: norm.cdf( (x - self.params.alpha - self.params.beta * self.params.distribution.ppf(u)) / np.abs(self.params.gamma) - ), - **params + ) ) - return rqmc.value, rqmc.error + return result.value, result.error - def compute_cdf(self, x: float, params: dict, integrator: Integrator = None) -> tuple[float, float]: + def compute_cdf(self, x: float, integrator: Integrator = RQMCIntegrator()) -> tuple[float, float]: """ Choose equation for cdf estimation depends on Mixture form Args: x (): point - params (): parameters of RQMC algorithm - integrator (): type of integrator to computing + integrator (): class of integrator with params to computing Returns: Computed pdf and error tolerance """ if isinstance(self.params, _NMMCanonicalDataCollector): - return self._canonical_compute_cdf(x, params, integrator) - return self._classical_compute_cdf(x, params, integrator) + return self._canonical_compute_cdf(x, integrator) + return self._classical_compute_cdf(x, integrator) - def _canonical_compute_pdf(self, x: float, params: dict, integrator: Integrator = None) -> tuple[float, float]: + def _canonical_compute_pdf(self, x: float, integrator: Integrator = RQMCIntegrator()) -> tuple[float, float]: """ Equation for canonical pdf Args: x (): point - params (): parameters of RQMC algorithm - integrator (): type of integrator to computing + integrator (): class of integrator with params to computing Returns: computed pdf and error tolerance """ - integrator = integrator or RQMCIntegrator() - rqmc = integrator.compute_integral(func= + result = integrator.compute(func= lambda u: (1 / np.abs(self.params.sigma)) * norm.pdf((x - self.params.distribution.ppf(u)) / np.abs(self.params.sigma)), - **params ) - return rqmc.value, rqmc.error + return result.value, result.error - def _classical_compute_pdf(self, x: float, params: dict, integrator: Integrator = None) -> tuple[float, float]: + def _classical_compute_pdf(self, x: float, integrator: Integrator = RQMCIntegrator()) -> tuple[float, float]: """ Equation for classic pdf Args: x (): point - params (): parameters of RQMC algorithm - integrator (): type of integrator to computing + integrator (): class of integrator with params to computing Returns: computed pdf and error tolerance """ - integrator = integrator or RQMCIntegrator() - rqmc = integrator.compute_integral(func= + result = integrator.compute(func= lambda u: (1 / np.abs(self.params.gamma)) * norm.pdf( (x - self.params.alpha - self.params.beta * self.params.distribution.ppf(u)) / np.abs(self.params.gamma) ), - **params ) - return rqmc.value, rqmc.error + return result.value, result.error - def compute_pdf(self, x: float, params: dict, integrator: Integrator = None) -> tuple[float, float]: + def compute_pdf(self, x: float, integrator: Integrator = RQMCIntegrator()) -> tuple[float, float]: """ Choose equation for pdf estimation depends on Mixture form Args: x (): point - params (): parameters of RQMC algorithm - integrator (): type of integrator to computing + integrator (): class of integrator with params to computing Returns: Computed pdf and error tolerance """ if isinstance(self.params, _NMMCanonicalDataCollector): - return self._canonical_compute_pdf(x, params, integrator) - return self._classical_compute_pdf(x, params, integrator) + return self._canonical_compute_pdf(x, integrator) + return self._classical_compute_pdf(x, integrator) def _classical_compute_log_pdf(self, x: float, params: dict) -> tuple[float, float]: """ diff --git a/src/mixtures/nmv_mixture.py b/src/mixtures/nmv_mixture.py index 9f2a7e0..a147cba 100644 --- a/src/mixtures/nmv_mixture.py +++ b/src/mixtures/nmv_mixture.py @@ -41,14 +41,13 @@ class NormalMeanVarianceMixtures(AbstractMixtures): def __init__(self, mixture_form: str, **kwargs: Any) -> None: super().__init__(mixture_form, **kwargs) - def _classical_moment(self, n: int, params: dict, integrator: Integrator = None) -> tuple[float, float]: + def _classical_moment(self, n: int, integrator: Integrator = RQMCIntegrator()) -> tuple[float, float]: """ Compute n-th moment of classical NMM Args: n (): Moment ordinal - params (): Parameters of integration algorithm - integrator (): type of integrator to computing + integrator (): class of integrator with params to computing Returns: moment approximation and error tolerance @@ -69,18 +68,16 @@ def integral_func(u: float) -> float: ) return result - integrator = integrator or RQMCIntegrator() - rqmc = integrator.compute_integral(func=lambda u: integral_func(u), **params) - return rqmc.value, rqmc.error + result = integrator.compute(func=lambda u: integral_func(u)) + return result.value, result.error - def _canonical_moment(self, n: int, params: dict, integrator: Integrator = None) -> tuple[float, float]: + def _canonical_moment(self, n: int, integrator: Integrator = RQMCIntegrator()) -> tuple[float, float]: """ Compute n-th moment of classical NMM Args: n (): Moment ordinal - params (): Parameters of integration algorithm - integrator (): type of integrator to computing + integrator (): class of integrator with params to computing Returns: moment approximation and error tolerance @@ -100,16 +97,15 @@ def integral_func(u: float) -> float: ) return result - integrator = integrator or RQMCIntegrator() - rqmc = integrator.compute_integral(func=lambda u: integral_func(u), **params) - return rqmc.value, rqmc.error + result = integrator.compute(func=lambda u: integral_func(u)) + return result.value, result.error - def compute_moment(self, n: int, params: dict, integrator: Integrator = None) -> tuple[float, float]: + def compute_moment(self, n: int, integrator: Integrator = RQMCIntegrator()) -> tuple[float, float]: if isinstance(self.params, _NMVMClassicDataCollector): - return self._classical_moment(n, params, integrator) - return self._canonical_moment(n, params, integrator) + return self._classical_moment(n, integrator) + return self._canonical_moment(n, integrator) - def _classical_cdf(self, x: float, params: dict, integrator: Integrator = None) -> tuple[float, float]: + def _classical_cdf(self, x: float, integrator: Integrator = RQMCIntegrator()) -> tuple[float, float]: def _inner_func(u: float) -> float: ppf = lru_cache()(self.params.distribution.ppf)(u) point = (x - self.params.alpha) / (np.sqrt(ppf) * self.params.gamma) - ( @@ -117,26 +113,24 @@ def _inner_func(u: float) -> float: ) return norm.cdf(point) - integrator = integrator or RQMCIntegrator() - rqmc = integrator.compute_integral(func=lambda u: _inner_func(u), **params) - return rqmc.value, rqmc.error + result = integrator.compute(func=lambda u: _inner_func(u)) + return result.value, result.error - def _canonical_cdf(self, x: float, params: dict, integrator: Integrator = None) -> tuple[float, float]: + def _canonical_cdf(self, x: float, integrator: Integrator = RQMCIntegrator()) -> tuple[float, float]: def _inner_func(u: float) -> float: ppf = self.params.distribution.ppf(u) point = (x - self.params.alpha) / (np.sqrt(ppf)) - (self.params.mu * np.sqrt(ppf)) return norm.cdf(point) - integrator = integrator or RQMCIntegrator() - rqmc = integrator.compute_integral(func=lambda u: _inner_func(u), **params) - return rqmc.value, rqmc.error + result = integrator.compute(func=lambda u: _inner_func(u)) + return result.value, result.error - def compute_cdf(self, x: float, params: dict, integrator: Integrator = None) -> tuple[float, float]: + def compute_cdf(self, x: float, integrator: Integrator = RQMCIntegrator()) -> tuple[float, float]: if isinstance(self.params, _NMVMClassicDataCollector): - return self._classical_cdf(x, params, integrator) - return self._canonical_cdf(x, params, integrator) + return self._classical_cdf(x, integrator) + return self._canonical_cdf(x, integrator) - def _classical_pdf(self, x: float, params: dict, integrator: Integrator = None) -> tuple[float, float]: + def _classical_pdf(self, x: float, integrator: Integrator = RQMCIntegrator()) -> tuple[float, float]: def _inner_func(u: float) -> float: ppf = self.params.distribution.ppf(u) return ( @@ -147,11 +141,10 @@ def _inner_func(u: float) -> float: ) ) - integrator = integrator or RQMCIntegrator() - rqmc = integrator.compute_integral(func=lambda u: _inner_func(u), **params) - return np.exp(self.params.beta * (x - self.params.alpha) / self.params.gamma**2) * rqmc.value, rqmc.error + result = integrator.compute(func=lambda u: _inner_func(u)) + return np.exp(self.params.beta * (x - self.params.alpha) / self.params.gamma**2) * result.value, result.error - def _canonical_pdf(self, x: float, params: dict, integrator: Integrator = None) -> tuple[float, float]: + def _canonical_pdf(self, x: float, integrator: Integrator = RQMCIntegrator()) -> tuple[float, float]: def _inner_func(u: float) -> float: ppf = self.params.distribution.ppf(u) return ( @@ -160,9 +153,8 @@ def _inner_func(u: float) -> float: * np.exp(-((x - self.params.alpha) ** 2 + self.params.mu**2 * ppf**2) / (2 * ppf)) ) - integrator = integrator or RQMCIntegrator() - rqmc = integrator.compute_integral(func=lambda u: _inner_func(u), **params) - return np.exp(self.params.mu * (x - self.params.alpha)) * rqmc.value, rqmc.error + result = integrator.compute(func=lambda u: _inner_func(u)) + return np.exp(self.params.mu * (x - self.params.alpha)) * result.value, result.error def _classical_log_pdf(self, x: float, params: dict) -> tuple[float, float]: def _inner_func(u: float) -> float: @@ -173,8 +165,8 @@ def _inner_func(u: float) -> float: + ppf * self.params.gamma**2 * np.log(2 * np.pi * ppf * self.params.gamma**2) ) / (2 * ppf * self.params.gamma**2) - rqmc = LogRQMC(lambda u: _inner_func(u), **params) - return rqmc() + result = LogRQMC(lambda u: _inner_func(u), **params) + return result() def _canonical_log_pdf(self, x: float, params: dict) -> tuple[float, float]: def _inner_func(u: float) -> float: @@ -183,13 +175,13 @@ def _inner_func(u: float) -> float: 2 * ppf ) - rqmc = LogRQMC(lambda u: _inner_func(u), **params) - return rqmc() + result = LogRQMC(lambda u: _inner_func(u), **params) + return result() - def compute_pdf(self, x: float, params: dict, integrator: Integrator = None) -> tuple[float, float]: + def compute_pdf(self, x: float, integrator: Integrator = RQMCIntegrator()) -> tuple[float, float]: if isinstance(self.params, _NMVMClassicDataCollector): - return self._classical_pdf(x, params, integrator) - return self._canonical_pdf(x, params, integrator) + return self._classical_pdf(x, integrator) + return self._canonical_pdf(x, integrator) def compute_logpdf(self, x: float, params: dict) -> tuple[Any, float]: if isinstance(self.params, _NMVMClassicDataCollector): diff --git a/src/mixtures/nv_mixture.py b/src/mixtures/nv_mixture.py index 1be3343..ffc4f08 100644 --- a/src/mixtures/nv_mixture.py +++ b/src/mixtures/nv_mixture.py @@ -40,13 +40,12 @@ class NormalVarianceMixtures(AbstractMixtures): def __init__(self, mixture_form: str, **kwargs: Any) -> None: super().__init__(mixture_form, **kwargs) - def compute_moment(self, n: int, params: dict, integrator: Integrator = None) -> tuple[float, float]: + def compute_moment(self, n: int, integrator: Integrator = RQMCIntegrator()) -> tuple[float, float]: """ Compute n-th moment of NVM Args: n (): Moment ordinal - params (): Parameters of integration algorithm - integrator (): type of integrator to computing + integrator (): class of integrator with params to computing Returns: moment approximation and error tolerance """ gamma = self.params.gamma if isinstance(self.params, _NVMClassicDataCollector) else 1 @@ -62,17 +61,15 @@ def integrate_func(u: float) -> float: for k in range(0, n + 1) ] ) - integrator = integrator or RQMCIntegrator() - result = integrator.compute_integral(func=integrate_func, **params) + result = integrator.compute(func=integrate_func) return result.value, result.error - def compute_cdf(self, x: float, params: dict, integrator: Integrator = None) -> tuple[float, float]: + def compute_cdf(self, x: float, integrator: Integrator = RQMCIntegrator()) -> tuple[float, float]: parametric_norm = norm(0, self.params.gamma if isinstance(self.params, _NVMClassicDataCollector) else 1) - integrator = integrator or RQMCIntegrator() - rqmc = integrator.compute_integral(func= - lambda u: parametric_norm.cdf((x - self.params.alpha) / np.sqrt(self.params.distribution.ppf(u))), **params + result = integrator.compute(func= + lambda u: parametric_norm.cdf((x - self.params.alpha) / np.sqrt(self.params.distribution.ppf(u))) ) - return rqmc.value, rqmc.error + return result.value, result.error @lru_cache() def _integrand_func(self, u: float, d: float, gamma: float) -> float: @@ -83,12 +80,11 @@ def _log_integrand_func(self, u: float, d: float, gamma: float | int | np.int64) ppf = self.params.distribution.ppf(u) return -(ppf * np.log(np.pi * 2 * ppf * gamma**2) + d) / (2 * ppf) - def compute_pdf(self, x: float, params: dict, integrator: Integrator = None) -> tuple[float, float]: + def compute_pdf(self, x: float, integrator: Integrator = RQMCIntegrator()) -> tuple[float, float]: gamma = self.params.gamma if isinstance(self.params, _NVMClassicDataCollector) else 1 d = (x - self.params.alpha) ** 2 / gamma**2 - integrator = integrator or RQMCIntegrator() - rqmc = integrator.compute_integral(func=lambda u: self._integrand_func(u, d, gamma), **params) - return rqmc.value, rqmc.error + result = integrator.compute(func=lambda u: self._integrand_func(u, d, gamma)) + return result.value, result.error def compute_logpdf(self, x: float, params: dict) -> tuple[float, float]: gamma = self.params.gamma if isinstance(self.params, _NVMClassicDataCollector) else 1