diff --git a/README.md b/README.md index 3b71d37..de58717 100644 --- a/README.md +++ b/README.md @@ -120,7 +120,7 @@ Computes given integral with the limits of integration being 0 and 1 using Rando ### Usage: ```Python -from src.algorithms.support_algorithms.rqmc import RQMC +from src.procedures.support.rqmc import RQMC rqmc = RQMC(lambda x: x**3 - x**2 + 1, error_tolerance=1e-5) ``` So `rqmc()[0]` is estimated integral value and `rqmc()[1]` is current error tolerance. diff --git a/jupiter_examples/nm_sigma_estimation_comparison.ipynb b/jupiter_examples/nm_sigma_estimation_comparison.ipynb index f986d13..10f695f 100644 --- a/jupiter_examples/nm_sigma_estimation_comparison.ipynb +++ b/jupiter_examples/nm_sigma_estimation_comparison.ipynb @@ -250,7 +250,7 @@ " \"\"\"\n", " generator = NMGenerator()\n", " mixture = NormalMeanMixtures(\"canonical\", sigma=real_sigma, distribution=distribution)\n", - " return generator.canonical_generate(mixture, sample_len)\n", + " return generator.generate(mixture, sample_len)\n", "\n", "def estimate_sigma_eigenvalue_based(sample, real_sigma, search_area, a, b):\n", " sample_len = len(sample)\n", diff --git a/src/algorithms/__init__.py b/src/algorithms/__init__.py deleted file mode 100644 index eaddba8..0000000 --- a/src/algorithms/__init__.py +++ /dev/null @@ -1,46 +0,0 @@ -from src.algorithms.semiparam_algorithms.nm_semi_param_algorithms.g_estimation_convolution import ( - NMSemiParametricGEstimation, -) -from src.algorithms.semiparam_algorithms.nm_semi_param_algorithms.sigma_estimation_eigenvalue_based import ( - SemiParametricMeanSigmaEstimationEigenvalueBased, -) -from src.algorithms.semiparam_algorithms.nm_semi_param_algorithms.sigma_estimation_empirical import ( - SemiParametricMeanSigmaEstimationEmpirical, -) -from src.algorithms.semiparam_algorithms.nv_semi_param_algorithms.g_estimation_given_mu import ( - SemiParametricNVEstimation, -) -from src.algorithms.semiparam_algorithms.nv_semi_param_algorithms.g_estimation_post_widder import ( - NVSemiParametricGEstimationPostWidder, -) -from src.algorithms.semiparam_algorithms.nvm_semi_param_algorithms.g_estimation_given_mu import ( - SemiParametricGEstimationGivenMu, -) -from src.algorithms.semiparam_algorithms.nvm_semi_param_algorithms.g_estimation_given_mu_rqmc_based import ( - SemiParametricGEstimationGivenMuRQMCBased, -) -from src.algorithms.semiparam_algorithms.nvm_semi_param_algorithms.g_estimation_post_widder import ( - SemiParametricGEstimationPostWidder, -) -from src.algorithms.semiparam_algorithms.nvm_semi_param_algorithms.mu_estimation import SemiParametricMuEstimation -from src.register.algorithm_purpose import AlgorithmPurpose -from src.register.register import Registry - -ALGORITHM_REGISTRY: Registry = Registry() -ALGORITHM_REGISTRY.register("mu_estimation", AlgorithmPurpose.NMV_SEMIPARAMETRIC)(SemiParametricMuEstimation) -ALGORITHM_REGISTRY.register("g_estimation_given_mu", AlgorithmPurpose.NMV_SEMIPARAMETRIC)( - SemiParametricGEstimationGivenMu -) -ALGORITHM_REGISTRY.register("g_estimation_convolution", AlgorithmPurpose.NM_SEMIPARAMETRIC)(NMSemiParametricGEstimation) -ALGORITHM_REGISTRY.register("sigma_estimation_eigenvalue_based", AlgorithmPurpose.NM_SEMIPARAMETRIC)( - SemiParametricMeanSigmaEstimationEigenvalueBased -) -ALGORITHM_REGISTRY.register("sigma_estimation_empirical", AlgorithmPurpose.NM_SEMIPARAMETRIC)( - SemiParametricMeanSigmaEstimationEmpirical -) -ALGORITHM_REGISTRY.register("g_estimation_given_mu_rqmc_based", AlgorithmPurpose.NMV_SEMIPARAMETRIC)( - SemiParametricGEstimationGivenMuRQMCBased -) -ALGORITHM_REGISTRY.register("g_estimation_post_widder", AlgorithmPurpose.NMV_SEMIPARAMETRIC)( - SemiParametricGEstimationPostWidder -) diff --git a/src/algorithms/support_algorithms/log_rqmc.py b/src/algorithms/support_algorithms/log_rqmc.py deleted file mode 100644 index 4423294..0000000 --- a/src/algorithms/support_algorithms/log_rqmc.py +++ /dev/null @@ -1,121 +0,0 @@ -from typing import Callable - -import numpy as np -import numpy._typing as tpg -import scipy - -from src.algorithms.support_algorithms.rqmc import RQMC - - -class LogRQMC(RQMC): - def __init__( - self, - func: Callable, - error_tolerance: float = 1e-6, - count: int = 25, - base_n: int = 2**6, - i_max: int = 100, - a: float = 0.00047, - ): - super().__init__(func, error_tolerance, count, base_n, i_max, a) - - @staticmethod - def lse(args: tpg.NDArray) -> float: - """ - Compute LSE - Args: - args (): Values - - Returns: LSE result - - """ - max_value = max(args) - return max_value + np.log(np.sum(np.exp(args - max_value))) - - def _independent_estimator(self, values: np._typing.NDArray) -> float: - """Apply function to row of matrix and find mean of row - - Args: - values: row of random values matrix - - Returns: mean of row - - """ - vfunc = np.vectorize(self.func) - return -np.log(len(values)) + self.lse(vfunc(values)) - - def _estimator(self, random_matrix: np._typing.NDArray) -> tuple[float, np._typing.NDArray]: - """Find mean of all rows - - Args: - random_matrix: matrix of random values - - Returns: mean of all rows - - """ - values = np.array(list(map(self._independent_estimator, random_matrix))) - return -np.log(self.count) + self.lse(values), values - - def _update_independent_estimator(self, i: int, old_value: float, new_values: np._typing.NDArray) -> float: - """Update mean of row - - Args: - i: step count - old_value: previous value of row on i-1 step - new_values: new generated row of random values - - Returns: Updated mean of row - - """ - - return -np.log(i + 1) + self.lse( - np.array(i * [old_value] + [self._independent_estimator(new_values[i * self.base_n :])]) - ) - - def _update( - self, j: int, old_values: np._typing.NDArray, random_matrix: np._typing.NDArray - ) -> tuple[float, np._typing.NDArray]: - """Update mean of all rows - - Args: - j: step count - old_values: previous values of row on i-1 step - random_matrix: new generated matrix of random values - - Returns:Updated mean of all rows - - """ - values = [] - for i in range(self.count): - old_value, new_values = old_values[i], random_matrix[i] - value = self._update_independent_estimator(j, old_value, new_values) - values.append(value) - np_values = np.array(values) - return -np.log(self.count) + self.lse(np_values), np_values - - def rqmc(self) -> tuple[float, float]: - """Main function of algorithm - - Returns: approximation for integral of function from 0 to 1 - - """ - sobol_sampler = scipy.stats.qmc.Sobol(d=1, scramble=False) - sobol_sample = np.repeat(sobol_sampler.random(self.base_n).transpose(), self.count, axis=0) - xor_sample = np.array(np.random.rand(1, self.count)[0]) - sample = self._digital_shift(sobol_sample, xor_sample) - approximation, values = self._estimator(sample) - current_error_tolerance = self._sigma(values, approximation) * self.z - for i in range(1, self.i_max): - if current_error_tolerance < self.error_tolerance: - return approximation, current_error_tolerance - sobol_sampler.reset() - sobol_sample = np.repeat(sobol_sampler.random(self.base_n * (i + 1)).transpose(), self.count, axis=0) - sample = self._digital_shift(sobol_sample, xor_sample) - approximation, values = self._update(i, values, sample) - current_error_tolerance = self._sigma(values, approximation) * self.z - return approximation, current_error_tolerance - - -if __name__ == "__main__": - log_rqmc = LogRQMC(lambda x: x**3 - x**2 + 1, i_max=100) - print(log_rqmc()) diff --git a/src/estimators/abstract_estimator.py b/src/estimators/abstract_estimator.py index 5afaea2..6bee8d3 100644 --- a/src/estimators/abstract_estimator.py +++ b/src/estimators/abstract_estimator.py @@ -2,7 +2,7 @@ from numpy import _typing -from src.algorithms import ALGORITHM_REGISTRY +from src.procedures import ALGORITHM_REGISTRY from src.estimators.estimate_result import EstimateResult from src.register.algorithm_purpose import AlgorithmPurpose @@ -14,7 +14,7 @@ class AbstractEstimator: algorithm_name: A string indicating chosen algorithm. params: A dictionary of algorithm parameters. estimate_result: Estimation result. - _registry: Registry that contains classes of all algorithms. + _registry: Registry that contains classes of all procedures. _purpose: Defines purpose of algorithm, one of the registry key. """ @@ -47,7 +47,7 @@ def set_params(self, algorithm_name: str, params: dict | None = None) -> Self: return self def get_available_algorithms(self) -> list[str]: - """Get all algorithms that can be used for current estimator class""" + """Get all procedures that can be used for current estimator class""" return [key[0] for key in self._registry.register_of_names.keys() if key[1] == self._purpose] def estimate(self, sample: _typing.ArrayLike) -> EstimateResult: diff --git a/src/estimators/semiparametric/abstract_semiparametric_estimator.py b/src/estimators/semiparametric/abstract_semiparametric_estimator.py index b26db4b..466f295 100644 --- a/src/estimators/semiparametric/abstract_semiparametric_estimator.py +++ b/src/estimators/semiparametric/abstract_semiparametric_estimator.py @@ -4,7 +4,7 @@ from src.estimators.estimate_result import EstimateResult -class AbstractSemiParametricEstimator(AbstractEstimator): +class AbstractSemiparEstim(AbstractEstimator): def __init__(self, algorithm_name: str, params: dict | None = None) -> None: super().__init__(algorithm_name, params) diff --git a/src/estimators/semiparametric/nm_semiparametric_estimator.py b/src/estimators/semiparametric/nm_semiparametric_estimator.py index 958ec5c..0513521 100644 --- a/src/estimators/semiparametric/nm_semiparametric_estimator.py +++ b/src/estimators/semiparametric/nm_semiparametric_estimator.py @@ -1,11 +1,11 @@ from numpy import _typing from src.estimators.estimate_result import EstimateResult -from src.estimators.semiparametric.abstract_semiparametric_estimator import AbstractSemiParametricEstimator +from src.estimators.semiparametric.abstract_semiparametric_estimator import AbstractSemiparEstim from src.register.algorithm_purpose import AlgorithmPurpose -class NMSemiParametricEstimator(AbstractSemiParametricEstimator): +class NMSemiparEstim(AbstractSemiparEstim): def __init__(self, algorithm_name: str, params: dict | None = None) -> None: super().__init__(algorithm_name, params) self._purpose = AlgorithmPurpose.NM_SEMIPARAMETRIC diff --git a/src/estimators/semiparametric/nmv_semiparametric_estimator.py b/src/estimators/semiparametric/nmv_semiparametric_estimator.py index e1f456d..93403e9 100644 --- a/src/estimators/semiparametric/nmv_semiparametric_estimator.py +++ b/src/estimators/semiparametric/nmv_semiparametric_estimator.py @@ -1,11 +1,11 @@ from numpy import _typing from src.estimators.estimate_result import EstimateResult -from src.estimators.semiparametric.abstract_semiparametric_estimator import AbstractSemiParametricEstimator +from src.estimators.semiparametric.abstract_semiparametric_estimator import AbstractSemiparEstim from src.register.algorithm_purpose import AlgorithmPurpose -class NMVSemiParametricEstimator(AbstractSemiParametricEstimator): +class NMVSemiparEstim(AbstractSemiparEstim): def __init__(self, algorithm_name: str, params: dict | None = None) -> None: super().__init__(algorithm_name, params) self._purpose = AlgorithmPurpose.NMV_SEMIPARAMETRIC diff --git a/src/estimators/semiparametric/nv_semiparametric_estimator.py b/src/estimators/semiparametric/nv_semiparametric_estimator.py index 95f0ce5..d2ec1d8 100644 --- a/src/estimators/semiparametric/nv_semiparametric_estimator.py +++ b/src/estimators/semiparametric/nv_semiparametric_estimator.py @@ -1,11 +1,11 @@ from numpy import _typing from src.estimators.estimate_result import EstimateResult -from src.estimators.semiparametric.abstract_semiparametric_estimator import AbstractSemiParametricEstimator +from src.estimators.semiparametric.abstract_semiparametric_estimator import AbstractSemiparEstim from src.register.algorithm_purpose import AlgorithmPurpose -class NVSemiParametricEstimator(AbstractSemiParametricEstimator): +class NVSemiparEstim(AbstractSemiparEstim): def __init__(self, algorithm_name: str, params: dict | None = None) -> None: super().__init__(algorithm_name, params) self._purpose = AlgorithmPurpose.NV_SEMIPARAMETRIC diff --git a/src/generators/nm_generator.py b/src/generators/nm_generator.py index 4c95b0f..5ecfe4a 100644 --- a/src/generators/nm_generator.py +++ b/src/generators/nm_generator.py @@ -9,7 +9,7 @@ class NMGenerator(AbstractGenerator): @staticmethod - def classical_generate(mixture: AbstractMixtures, size: int) -> tpg.NDArray: + def generate(mixture: AbstractMixtures, size: int) -> tpg.NDArray: """Generate a sample of given size. Classical form of NMM Args: @@ -27,25 +27,6 @@ def classical_generate(mixture: AbstractMixtures, size: int) -> tpg.NDArray: raise ValueError("Mixture must be NormalMeanMixtures") mixing_values = mixture.params.distribution.rvs(size=size) normal_values = scipy.stats.norm.rvs(size=size) + if mixture.mixture_form == "canonical": + return mixing_values + mixture.params.sigma * normal_values return mixture.params.alpha + mixture.params.beta * mixing_values + mixture.params.gamma * normal_values - - @staticmethod - def canonical_generate(mixture: AbstractMixtures, size: int) -> tpg.NDArray: - """Generate a sample of given size. Canonical form of NMM - - Args: - mixture: Normal Mean Mixture - size: length of sample - - Returns: sample of given size - - Raises: - ValueError: If mixture is not a Normal Mean Mixture - - """ - - if not isinstance(mixture, NormalMeanMixtures): - raise ValueError("Mixture must be NormalMeanMixtures") - mixing_values = mixture.params.distribution.rvs(size=size) - normal_values = scipy.stats.norm.rvs(size=size) - return mixing_values + mixture.params.sigma * normal_values diff --git a/src/generators/nmv_generator.py b/src/generators/nmv_generator.py index 5ff3221..8270a1a 100644 --- a/src/generators/nmv_generator.py +++ b/src/generators/nmv_generator.py @@ -9,7 +9,7 @@ class NMVGenerator(AbstractGenerator): @staticmethod - def classical_generate(mixture: AbstractMixtures, size: int) -> tpg.NDArray: + def generate(mixture: AbstractMixtures, size: int) -> tpg.NDArray: """Generate a sample of given size. Classical form of NMVM Args: @@ -27,29 +27,10 @@ def classical_generate(mixture: AbstractMixtures, size: int) -> tpg.NDArray: raise ValueError("Mixture must be NormalMeanMixtures") mixing_values = mixture.params.distribution.rvs(size=size) normal_values = scipy.stats.norm.rvs(size=size) + if mixture.mixture_form == "canonical": + return mixture.params.alpha + mixture.params.mu * mixing_values + (mixing_values ** 0.5) * normal_values return ( mixture.params.alpha + mixture.params.beta * mixing_values + mixture.params.gamma * (mixing_values**0.5) * normal_values - ) - - @staticmethod - def canonical_generate(mixture: AbstractMixtures, size: int) -> tpg.NDArray: - """Generate a sample of given size. Canonical form of NMVM - - Args: - mixture: Normal Mean Variance Mixtures - size: length of sample - - Returns: sample of given size - - Raises: - ValueError: If mixture type is not Normal Mean Variance Mixtures - - """ - - if not isinstance(mixture, NormalMeanVarianceMixtures): - raise ValueError("Mixture must be NormalMeanMixtures") - mixing_values = mixture.params.distribution.rvs(size=size) - normal_values = scipy.stats.norm.rvs(size=size) - return mixture.params.alpha + mixture.params.mu * mixing_values + (mixing_values**0.5) * normal_values + ) \ No newline at end of file diff --git a/src/generators/nv_generator.py b/src/generators/nv_generator.py index faa55f5..712e918 100644 --- a/src/generators/nv_generator.py +++ b/src/generators/nv_generator.py @@ -9,7 +9,7 @@ class NVGenerator(AbstractGenerator): @staticmethod - def classical_generate(mixture: AbstractMixtures, size: int) -> tpg.NDArray: + def generate(mixture: AbstractMixtures, size: int) -> tpg.NDArray: """Generate a sample of given size. Classical form of NVM Args: @@ -27,25 +27,6 @@ def classical_generate(mixture: AbstractMixtures, size: int) -> tpg.NDArray: raise ValueError("Mixture must be NormalMeanMixtures") mixing_values = mixture.params.distribution.rvs(size=size) normal_values = scipy.stats.norm.rvs(size=size) - return mixture.params.alpha + mixture.params.gamma * (mixing_values**0.5) * normal_values - - @staticmethod - def canonical_generate(mixture: AbstractMixtures, size: int) -> tpg.NDArray: - """Generate a sample of given size. Canonical form of NVM - - Args: - mixture: Normal Variance Mixtures - size: length of sample - - Returns: sample of given size - - Raises: - ValueError: If mixture type is not Normal Variance Mixtures - - """ - - if not isinstance(mixture, NormalVarianceMixtures): - raise ValueError("Mixture must be NormalMeanMixtures") - mixing_values = mixture.params.distribution.rvs(size=size) - normal_values = scipy.stats.norm.rvs(size=size) - return mixture.params.alpha + (mixing_values**0.5) * normal_values + if mixture.mixture_form == "canonical": + return mixture.params.alpha + (mixing_values ** 0.5) * normal_values + return mixture.params.alpha + mixture.params.gamma * (mixing_values**0.5) * normal_values \ No newline at end of file diff --git a/src/mixtures/abstract_mixture.py b/src/mixtures/abstract_mixture.py index 13c8e42..de31713 100644 --- a/src/mixtures/abstract_mixture.py +++ b/src/mixtures/abstract_mixture.py @@ -1,10 +1,13 @@ from abc import ABCMeta, abstractmethod from dataclasses import fields -from typing import Any +from typing import Any, List, Tuple, Union, Dict, Type +import numpy as np +from numpy.typing import NDArray from scipy.stats import rv_continuous from scipy.stats.distributions import rv_frozen +from src.procedures.support.integrator import Integrator class AbstractMixtures(metaclass=ABCMeta): """Base class for Mixtures""" @@ -12,14 +15,21 @@ class AbstractMixtures(metaclass=ABCMeta): _classical_collector: Any _canonical_collector: Any - @abstractmethod - def __init__(self, mixture_form: str, **kwargs: Any) -> None: + def __init__( + self, + mixture_form: str, + **kwargs: Any + ) -> None: """ - Args: - mixture_form: Form of Mixture classical or Canonical - **kwargs: Parameters of Mixture + mixture_form: Form of Mixture classical or canonical + integrator_cls: Class implementing Integrator protocol (default: RQMCIntegrator) + integrator_params: Parameters for integrator constructor (default: {{}}) + **kwargs: Parameters of Mixture (alpha, gamma, etc.) """ + self.mixture_form = mixture_form + + if mixture_form == "classical": self.params = self._params_validation(self._classical_collector, kwargs) elif mixture_form == "canonical": @@ -28,40 +38,92 @@ 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]: + ... + + def compute_moment( + self, + x: Union[List[int], int, NDArray[np.float64]], + integrator: Integrator + ) -> Union[List[Tuple[float, float]], Tuple[float, float], NDArray[Any]]: + if isinstance(x, np.ndarray): + return np.array([self._compute_moment(xp, integrator) for xp in x], dtype=object) + elif isinstance(x, list): + return [self._compute_moment(xp, integrator) for xp in x] + elif isinstance(x, int): + return self._compute_moment(x, integrator) + else: + raise TypeError(f"Unsupported type for x: {type(x)}") @abstractmethod - def compute_cdf(self, x: float, params: dict) -> tuple[float, float]: ... + def _compute_pdf(self, x: float, integrator: Integrator) -> Tuple[float, float]: + ... + + def compute_pdf( + self, + x: Union[List[float], float, NDArray[np.float64]], + integrator: Integrator + ) -> Union[List[Tuple[float, float]], Tuple[float, float], NDArray[Any]]: + if isinstance(x, np.ndarray): + return np.array([self._compute_pdf(xp, integrator) for xp in x], dtype=object) + elif isinstance(x, list): + return [self._compute_pdf(xp, integrator) for xp in x] + elif isinstance(x, float): + return self._compute_pdf(x, integrator) + else: + raise TypeError(f"Unsupported type for x: {type(x)}") @abstractmethod - def compute_pdf(self, x: float, params: dict) -> tuple[float, float]: ... + def _compute_logpdf(self, x: float, integrator: Integrator) -> Tuple[float, float]: + ... + + def compute_logpdf( + self, + x: Union[List[float], float, NDArray[np.float64]], + integrator: Integrator + ) -> Union[List[Tuple[float, float]], Tuple[float, float], NDArray[Any]]: + if isinstance(x, np.ndarray): + return np.array([self._compute_logpdf(xp, integrator) for xp in x], dtype=object) + elif isinstance(x, list): + return [self._compute_logpdf(xp, integrator) for xp in x] + elif isinstance(x, float): + return self._compute_logpdf(x, integrator) + else: + raise TypeError(f"Unsupported type for x: {type(x)}") @abstractmethod - def compute_logpdf(self, x: float, params: dict) -> tuple[float, float]: ... - - def _params_validation(self, data_collector: Any, params: dict[str, float | rv_continuous | rv_frozen]) -> Any: - """Mixture Parameters Validation - - Args: - data_collector: Dataclass that collect parameters of Mixture - params: Input parameters - - Returns: Instance of dataclass - - Raises: - ValueError: If given parameters is unexpected - ValueError: If parameter type is invalid - ValueError: If parameters age not given - - """ - + def _compute_cdf(self, x: float, integrator: Integrator) -> Tuple[float, float]: + ... + + def compute_cdf( + self, + x: Union[List[float], float, NDArray[np.float64]], + integrator: Integrator + ) -> Union[List[Tuple[float, float]], Tuple[float, float], NDArray[Any]]: + if isinstance(x, np.ndarray): + return np.array([self._compute_cdf(xp, integrator) for xp in x], dtype=object) + elif isinstance(x, list): + return [self._compute_cdf(xp, integrator) for xp in x] + elif isinstance(x, float): + return self._compute_cdf(x, integrator) + else: + raise TypeError(f"Unsupported type for x: {type(x)}") + + def _params_validation( + self, + data_collector: Any, + params: dict[str, float | rv_continuous | rv_frozen] + ) -> Any: + """Mixture Parameters Validation""" dataclass_fields = fields(data_collector) if len(params) != len(dataclass_fields): raise ValueError(f"Expected {len(dataclass_fields)} arguments, got {len(params)}") - names_and_types = dict((field.name, field.type) for field in dataclass_fields) - for pair in params.items(): - if pair[0] not in names_and_types: - raise ValueError(f"Unexpected parameter {pair[0]}") - if not isinstance(pair[1], names_and_types[pair[0]]): - raise ValueError(f"Type missmatch: {pair[0]} should be {names_and_types[pair[0]]}, not {type(pair[1])}") + names_and_types = {field.name: field.type for field in dataclass_fields} + for name, value in params.items(): + if name not in names_and_types: + raise ValueError(f"Unexpected parameter {name}") + if not isinstance(value, names_and_types[name]): + raise ValueError( + f"Type mismatch: {name} should be {names_and_types[name]}, not {type(value)}" + ) return data_collector(**params) diff --git a/src/mixtures/nm_mixture.py b/src/mixtures/nm_mixture.py index 33c7580..1f1e15f 100644 --- a/src/mixtures/nm_mixture.py +++ b/src/mixtures/nm_mixture.py @@ -1,261 +1,106 @@ from dataclasses import dataclass -from typing import Any +from typing import Any, Type, Dict, Tuple 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.procedures.support.integrator import Integrator +from src.procedures.support.rqmc import RQMCIntegrator +from src.procedures.support.rqmc import RQMC +from src.procedures.support.log_rqmc import LogRQMCIntegrator +from src.procedures.support.quad_integrator import QuadIntegrator from src.mixtures.abstract_mixture import AbstractMixtures - @dataclass class _NMMClassicDataCollector: - """TODO: Change typing from float | int | etc to Protocol with __addition__ __multiplication__ __subtraction__""" - - """Data Collector for parameters of classical NMM""" alpha: float | int | np.int64 - beta: float | int | np.int64 + beta: float | int | np.int64 gamma: float | int | np.int64 distribution: rv_frozen | rv_continuous - @dataclass class _NMMCanonicalDataCollector: - """TODO: Change typing from float | int | etc to Protocol with __addition__ __multiplication__ __subtraction__""" - - """Data Collector for parameters of canonical NMM""" sigma: float | int | np.int64 distribution: rv_frozen | rv_continuous - class NormalMeanMixtures(AbstractMixtures): _classical_collector = _NMMClassicDataCollector _canonical_collector = _NMMCanonicalDataCollector - def __init__(self, mixture_form: str, **kwargs: Any) -> None: - """ - Read Doc of Parent Method - """ - + def __init__( + self, + mixture_form: str, + **kwargs: Any + ) -> None: super().__init__(mixture_form, **kwargs) def _params_validation(self, data_collector: Any, params: dict[str, float | rv_continuous | rv_frozen]) -> Any: - """ - Read parent method doc - - Raises: - ValueError: If canonical Mixture has negative sigma parameter - - """ - data_class = super()._params_validation(data_collector, params) if hasattr(data_class, "sigma") and data_class.sigma <= 0: - raise ValueError("Sigma cant be zero or negative") + raise ValueError("Sigma can't be zero or negative") if hasattr(data_class, "gamma") and data_class.gamma == 0: - raise ValueError("Gamma cant be zero") + raise ValueError("Gamma can't be zero") return data_class - def _classical_moment(self, n: int, params: dict) -> tuple[float, float]: - """ - Compute n-th moment of classical NMM - - Args: - n (): Moment ordinal - params (): Parameters of integration algorithm - - Returns: moment approximation and error tolerance - - """ - mixture_moment = 0 - error_tolerance = 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 = 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 - - def _canonical_moment(self, n: int, params: dict) -> tuple[float, float]: - """ - Compute n-th moment of canonical NMM - - Args: - n (): Moment ordinal - params (): Parameters of integration algorithm - - Returns: moment approximation and error tolerance - - """ - mixture_moment = 0 - error_tolerance = 0 - 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) - return mixture_moment, error_tolerance - - def compute_moment(self, n: int, params: dict) -> tuple[float, float]: - """ - Compute n-th moment of NMM - - Args: - n (): Moment ordinal - params (): Parameters of integration algorithm - - Returns: moment approximation and error tolerance - - """ - if isinstance(self.params, _NMMClassicDataCollector): - return self._classical_moment(n, params) - return self._canonical_moment(n, params) - - def _canonical_compute_cdf(self, x: float, params: dict) -> tuple[float, float]: - """ - Equation for canonical cdf - Args: - x (): point - params (): parameters of RQMC algorithm - - 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() - - def _classical_compute_cdf(self, x: float, params: dict) -> tuple[float, float]: - """ - Equation for classic cdf - Args: - x (): point - params (): parameters of RQMC algorithm - - Returns: computed cdf and error tolerance - - """ - rqmc = RQMC( - lambda u: norm.cdf( - (x - self.params.alpha - self.params.beta * self.params.distribution.ppf(u)) / np.abs(self.params.gamma) - ), - **params - ) - return rqmc() - - def compute_cdf(self, x: float, params: dict) -> tuple[float, float]: - """ - Choose equation for cdf estimation depends on Mixture form - Args: - x (): point - params (): parameters of RQMC algorithm - - 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) - - def _canonical_compute_pdf(self, x: float, params: dict) -> tuple[float, float]: - """ - Equation for canonical pdf - Args: - x (): point - params (): parameters of RQMC algorithm - - Returns: computed pdf and error tolerance - - """ - rqmc = RQMC( - lambda u: (1 / np.abs(self.params.sigma)) - * norm.pdf((x - self.params.distribution.ppf(u)) / np.abs(self.params.sigma)), - **params - ) - return rqmc() - - def _classical_compute_pdf(self, x: float, params: dict) -> tuple[float, float]: - """ - Equation for classic pdf - Args: - x (): point - params (): parameters of RQMC algorithm - - Returns: computed pdf and error tolerance - - """ - rqmc = RQMC( - 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() - - def compute_pdf(self, x: float, params: dict) -> tuple[float, float]: - """ - Choose equation for pdf estimation depends on Mixture form - Args: - x (): point - params (): parameters of RQMC algorithm - - 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) - - def _classical_compute_log_pdf(self, x: float, params: dict) -> tuple[float, float]: - """ - Equation for classic log pdf - Args: - x (): point - params (): parameters of LogRQMC algorithm - - Returns: computed log pdf and error tolerance - - """ - rqmc = LogRQMC( - lambda u: np.log(1 / np.abs(self.params.gamma)) - + norm.logpdf( - (x - self.params.alpha - self.params.beta * self.params.distribution.ppf(u)) / np.abs(self.params.gamma) - ), - **params - ) - return rqmc() - - def _canonical_compute_log_pdf(self, x: float, params: dict) -> tuple[float, float]: - """ - Equation for canonical log pdf - Args: - x (): point - params (): parameters of LogRQMC algorithm - - Returns: computed log pdf and error tolerance - - """ - rqmc = LogRQMC( - lambda u: np.log(1 / np.abs(self.params.sigma)) - + norm.logpdf((x - self.params.distribution.ppf(u)) / np.abs(self.params.sigma)), - **params - ) - return rqmc() - - def compute_logpdf(self, x: float, params: dict) -> tuple[float, float]: - """ - Choose equation for log pdf estimation depends on Mixture form - Args: - x (): point - params (): parameters of LogRQMC algorithm - - Returns: Computed log pdf and error tolerance - - """ - if isinstance(self.params, _NMMCanonicalDataCollector): - return self._canonical_compute_log_pdf(x, params) - return self._classical_compute_log_pdf(x, params) + def _compute_moment(self, n: int, integrator: Integrator=QuadIntegrator()) -> Tuple[float, float]: + mixture_moment = 0.0 + error = 0.0 + if self.mixture_form == "classical": + for k in range(n + 1): + for l in range(k + 1): + coeff = binom(n, n - k) * binom(k, k - l) + def mix(u: float) -> float: + return ( + self.params.distribution.ppf(u) ** (k - l) + ) + + res = integrator.compute(mix) + mixture_moment += coeff * (self.params.beta ** (k - l)) * (self.params.gamma ** l) * (self.params.alpha ** (n - k)) * res.value * norm.moment(l) + error += coeff * (self.params.beta ** (k - l)) * (self.params.gamma ** l) * (self.params.alpha ** (n - k)) * res.error * norm.moment(l) + else: + for k in range(n + 1): + coeff = binom(n, n - k) + def mix(u: float) -> float: + return self.params.distribution.ppf(u) ** (n - k) + res = integrator.compute(mix) + mixture_moment += coeff * (self.params.sigma ** k) * res.value * norm.moment(k) + error += coeff * (self.params.sigma ** k) * res.error * norm.moment(k) + return mixture_moment, error + + def _compute_cdf(self, x: float, integrator: Integrator=RQMCIntegrator()) -> Tuple[float, float]: + if self.mixture_form == "classical": + def fn(u: float) -> float: + p = self.params.distribution.ppf(u) + return norm.cdf((x - self.params.alpha - self.params.beta * p) / abs(self.params.gamma)) + else: + def fn(u: float) -> float: + p = self.params.distribution.ppf(u) + return norm.cdf((x - p) / abs(self.params.sigma)) + res = integrator.compute(fn) + return res.value, res.error + + def _compute_pdf(self, x: float, integrator: Integrator=RQMCIntegrator()) -> Tuple[float, float]: + if self.mixture_form == "classical": + def fn(u: float) -> float: + p = self.params.distribution.ppf(u) + return (1 / abs(self.params.gamma)) * norm.pdf((x - self.params.alpha - self.params.beta * p) / abs(self.params.gamma)) + else: + def fn(u: float) -> float: + p = self.params.distribution.ppf(u) + return (1 / abs(self.params.sigma)) * norm.pdf((x - p) / abs(self.params.sigma)) + res = integrator.compute(fn) + return res.value, res.error + + def _compute_logpdf(self, x: float, integrator: Integrator=LogRQMCIntegrator()) -> Tuple[float, float]: + if self.mixture_form == "classical": + def fn(u: float) -> float: + p = self.params.distribution.ppf(u) + return np.log(1 / abs(self.params.gamma)) + norm.logpdf((x - self.params.alpha - self.params.beta * p) / abs(self.params.gamma)) + else: + def fn(u: float) -> float: + p = self.params.distribution.ppf(u) + return np.log(1 / abs(self.params.sigma)) + norm.logpdf((x - p) / abs(self.params.sigma)) + res = integrator.compute(fn) + return res.value, res.error diff --git a/src/mixtures/nmv_mixture.py b/src/mixtures/nmv_mixture.py index c090ff3..6bd4dcb 100644 --- a/src/mixtures/nmv_mixture.py +++ b/src/mixtures/nmv_mixture.py @@ -1,191 +1,112 @@ from dataclasses import dataclass from functools import lru_cache -from typing import Any +from typing import Any, Type, Dict, Tuple import numpy as np from scipy.special import binom -from scipy.stats import geninvgauss, norm, rv_continuous +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.procedures.support.integrator import Integrator +from src.procedures.support.rqmc import RQMCIntegrator +from src.procedures.support.rqmc import RQMC +from src.procedures.support.log_rqmc import LogRQMCIntegrator from src.mixtures.abstract_mixture import AbstractMixtures - @dataclass class _NMVMClassicDataCollector: - """TODO: Change typing from float | int | etc to Protocol with __addition__ __multiplication__ __subtraction__""" - - """Data Collector for parameters of classical NMVM""" alpha: float | int | np.int64 - beta: float | int | np.int64 + beta: float | int | np.int64 gamma: float | int | np.int64 distribution: rv_frozen | rv_continuous - @dataclass class _NMVMCanonicalDataCollector: - """TODO: Change typing from float | int | etc to Protocol with __addition__ __multiplication__ __subtraction__""" - - """Data Collector for parameters of canonical NMVM""" alpha: float | int | np.int64 - mu: float | int | np.int64 + mu: float | int | np.int64 distribution: rv_frozen | rv_continuous - class NormalMeanVarianceMixtures(AbstractMixtures): _classical_collector = _NMVMClassicDataCollector _canonical_collector = _NMVMCanonicalDataCollector - def __init__(self, mixture_form: str, **kwargs: Any) -> None: + 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]: - """ - Compute n-th moment of classical NMM - - Args: - n (): Moment ordinal - params (): Parameters of integration algorithm - - Returns: moment approximation and error tolerance - - """ - - def integral_func(u: float) -> float: - result = 0 - for k in range(0, n + 1): - for l in range(0, k + 1): - result += ( - binom(n, n - k) - * binom(k, k - l) - * (self.params.beta ** (k - l)) - * (self.params.gamma**l) - * self.params.distribution.ppf(u) ** (k - l / 2) - * (self.params.alpha ** (n - k)) - * norm.moment(l) - ) - return result - - rqmc = RQMC(lambda u: integral_func(u), **params) - return rqmc() - - def _canonical_moment(self, n: int, params: dict) -> tuple[float, float]: - """ - Compute n-th moment of classical NMM - - Args: - n (): Moment ordinal - params (): Parameters of integration algorithm - - Returns: moment approximation and error tolerance - - """ - - def integral_func(u: float) -> float: - result = 0 - for k in range(0, n + 1): - for l in range(0, k + 1): - result += ( - binom(n, n - k) - * binom(k, k - l) - * (self.params.nu ** (k - l)) - * self.params.distribution.ppf(u) ** (k - l / 2) - * (self.params.alpha ** (n - k)) - * norm.moment(l) - ) - return result - - rqmc = RQMC(lambda u: integral_func(u), **params) - return rqmc() - - def compute_moment(self, n: int, params: dict) -> tuple[float, float]: - if isinstance(self.params, _NMVMClassicDataCollector): - return self._classical_moment(n, params) - return self._canonical_moment(n, params) - - def _classical_cdf(self, x: float, params: dict) -> 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) - ( - self.params.beta / self.params.gamma * np.sqrt(ppf) - ) - return norm.cdf(point) - - rqmc = RQMC(lambda u: _inner_func(u), **params) - return rqmc() - - def _canonical_cdf(self, x: float, params: dict) -> 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() - - def compute_cdf(self, x: float, params: dict) -> tuple[float, float]: - if isinstance(self.params, _NMVMClassicDataCollector): - return self._classical_cdf(x, params) - return self._canonical_cdf(x, params) - - def _classical_pdf(self, x: float, params: dict) -> tuple[float, float]: - def _inner_func(u: float) -> float: - ppf = self.params.distribution.ppf(u) - return ( - 1 - / np.sqrt(2 * np.pi * ppf * self.params.gamma**2) - * np.exp( - -((x - self.params.alpha) ** 2 + self.params.beta**2 * ppf**2) / (2 * ppf * self.params.gamma**2) + def _compute_moment(self, n: int, integrator: Integrator=RQMCIntegrator()) -> Tuple[float, float]: + gamma = getattr(self.params, 'gamma', None) + + def integrand(u: float) -> float: + s = 0.0 + for k in range(n + 1): + for l in range(k + 1): + if self.mixture_form == 'classical': + term = ( + binom(n, n - k) + * binom(k, k - l) + * (self.params.beta ** (k - l)) + * (self.params.gamma ** l) + * (self.params.distribution.ppf(u) ** (k - l/2)) + * norm.moment(l) + ) + else: + term = ( + binom(n, n - k) + * binom(k, k - l) + * (self.params.mu ** (k - l)) + * (self.params.distribution.ppf(u) ** (k - l/2)) + * norm.moment(l) + ) + s += term + return s + + res = integrator.compute(integrand) + return res.value, res.error + + def _compute_cdf(self, x: float, integrator: Integrator=RQMCIntegrator()) -> Tuple[float, float]: + def integrand(u: float) -> float: + p = self.params.distribution.ppf(u) + if self.mixture_form == 'classical': + return norm.cdf((x - self.params.alpha) / (np.sqrt(p) * self.params.gamma)) + return norm.cdf((x - self.params.alpha) / np.sqrt(p) - self.params.mu * np.sqrt(p)) + + res = integrator.compute(integrand) + return res.value, res.error + + def _compute_pdf(self, x: float, integrator: Integrator=RQMCIntegrator()) -> Tuple[float, float]: + def integrand(u: float) -> float: + p = self.params.distribution.ppf(u) + if self.mixture_form == 'classical': + return ( + 1 / np.sqrt(2 * np.pi * p * self.params.gamma ** 2) + * np.exp(-((x - self.params.alpha) ** 2 + self.params.beta ** 2 * p ** 2) / (2 * p * self.params.gamma ** 2)) ) - ) - - 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] - - def _canonical_pdf(self, x: float, params: dict) -> tuple[float, float]: - def _inner_func(u: float) -> float: - ppf = self.params.distribution.ppf(u) return ( - 1 - / np.sqrt(2 * np.pi * ppf) - * 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] - - def _classical_log_pdf(self, x: float, params: dict) -> tuple[float, float]: - def _inner_func(u: float) -> float: - ppf = self.params.distribution.ppf(u) - return -( - (x - self.params.alpha) ** 2 - + ppf**2 * self.params.beta**2 - + 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() - - def _canonical_log_pdf(self, x: float, params: dict) -> tuple[float, float]: - def _inner_func(u: float) -> float: - ppf = self.params.distribution.ppf(u) - return -((x - self.params.alpha) ** 2 + ppf**2 * self.params.mu**2 + ppf * np.log(2 * np.pi * ppf)) / ( - 2 * ppf + 1 / np.sqrt(2 * np.pi * p) + * np.exp(-((x - self.params.alpha) ** 2 + self.params.mu ** 2 * p ** 2) / (2 * p)) ) - rqmc = LogRQMC(lambda u: _inner_func(u), **params) - return rqmc() - - def compute_pdf(self, x: float, params: dict) -> tuple[float, float]: - if isinstance(self.params, _NMVMClassicDataCollector): - return self._classical_pdf(x, params) - return self._canonical_pdf(x, params) - - def compute_logpdf(self, x: float, params: dict) -> tuple[Any, float]: - if isinstance(self.params, _NMVMClassicDataCollector): - int_result = self._classical_log_pdf(x, params) - return self.params.beta * (x - self.params.alpha) / self.params.gamma**2 + int_result[0], int_result[1] - int_result = self._canonical_log_pdf(x, params) - return self.params.mu * (x - self.params.alpha) + int_result[0], int_result[1] + res = integrator.compute(integrand) + if self.mixture_form == 'classical': + val = np.exp(self.params.beta * (x - self.params.alpha) / self.params.gamma ** 2) * res.value + else: + val = np.exp(self.params.mu * (x - self.params.alpha)) * res.value + return val, res.error + + def _compute_logpdf(self, x: float, integrator: Integrator=LogRQMCIntegrator()) -> Tuple[float, float]: + def integrand(u: float) -> float: + p = self.params.distribution.ppf(u) + if self.mixture_form == 'classical': + return -((x - self.params.alpha) ** 2 + p ** 2 * self.params.beta ** 2 + p * self.params.gamma ** 2 * np.log(2 * np.pi * p * self.params.gamma ** 2)) / (2 * p * self.params.gamma ** 2) + return -((x - self.params.alpha) ** 2 + p ** 2 * self.params.mu ** 2 + p * np.log(2 * np.pi * p)) / (2 * p) + + res = integrator.compute(integrand) + if self.mixture_form == 'classical': + val = self.params.beta * (x - self.params.alpha) / self.params.gamma ** 2 + res.value + else: + val = self.params.mu * (x - self.params.alpha) + res.value + return val, res.error diff --git a/src/mixtures/nv_mixture.py b/src/mixtures/nv_mixture.py index 7905e6e..20f75ae 100644 --- a/src/mixtures/nv_mixture.py +++ b/src/mixtures/nv_mixture.py @@ -1,93 +1,97 @@ from dataclasses import dataclass from functools import lru_cache -from typing import Any +from typing import Any, Type, Dict import numpy as np 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.mixtures.abstract_mixture import AbstractMixtures +from src.procedures.support.log_rqmc import LogRQMCIntegrator +from src.procedures.support.integrator import Integrator +from src.procedures.support.rqmc import RQMCIntegrator +from src.procedures.support.rqmc import RQMC +from src.procedures.support.quad_integrator import QuadIntegrator +from src.mixtures.abstract_mixture import AbstractMixtures @dataclass class _NVMClassicDataCollector: - """TODO: Change typing from float | int | etc to Protocol with __addition__ __multiplication__ __subtraction__""" - - """Data Collector for parameters of classical NVM""" alpha: float | int | np.int64 gamma: float | int | np.int64 distribution: rv_frozen | rv_continuous - @dataclass class _NVMCanonicalDataCollector: - """TODO: Change typing from float | int | etc to Protocol with __addition__ __multiplication__ __subtraction__""" - - """Data Collector for parameters of canonical NVM""" alpha: float | int | np.int64 distribution: rv_frozen | rv_continuous - class NormalVarianceMixtures(AbstractMixtures): - _classical_collector = _NVMClassicDataCollector _canonical_collector = _NVMCanonicalDataCollector - 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]: - """ - Compute n-th moment of NVM - Args: - n (): Moment ordinal - params (): Parameters of integration algorithm - Returns: moment approximation and error tolerance - """ - gamma = self.params.gamma if isinstance(self.params, _NVMClassicDataCollector) else 1 - - def integrate_func(u: float) -> float: + def __init__( + self, + mixture_form: str, + integrator_cls: Type[Integrator] = RQMCIntegrator, + integrator_params: Dict[str, Any] = None, + **kwargs: Any + ) -> None: + super().__init__(mixture_form, integrator_cls=integrator_cls, integrator_params=integrator_params, **kwargs) + self.integrator_cls = integrator_cls + self.integrator_params = integrator_params or {} + + def _compute_moment(self, n: int, integrator: Integrator=QuadIntegrator()) -> tuple[float, float]: + gamma = getattr(self.params, 'gamma', 1) + + def integrand(u: float) -> float: return sum( - [ - binom(n, k) - * (gamma**k) - * (self.params.alpha ** (n - k)) - * (self.params.distribution.ppf(u) ** (k / 2)) - * norm.moment(k) - for k in range(0, n + 1) - ] + binom(n, k) + * (gamma ** k) + * (self.params.alpha ** (n - k)) + * (self.params.distribution.ppf(u) ** (k / 2)) + * norm.moment(k) + for k in range(n + 1) ) - result = RQMC(integrate_func, **params)() - return result + result = integrator.compute(integrand) + return result.value, result.error + + def _compute_cdf(self, x: float, integrator: Integrator=QuadIntegrator()) -> tuple[float, float]: + gamma = getattr(self.params, 'gamma', 1) + param_norm = norm(0, gamma) + + def integrand(u: float) -> float: + return param_norm.cdf((x - self.params.alpha) / np.sqrt(self.params.distribution.ppf(u))) + + result = integrator.compute(integrand) + return result.value, result.error + + def _compute_pdf(self, x: float, integrator: Integrator=QuadIntegrator()) -> tuple[float, float]: + gamma = getattr(self.params, 'gamma', 1) + d = (x - self.params.alpha) ** 2 / gamma ** 2 + + def integrand(u: float) -> float: + return self._integrand_func(u, d, gamma) + + result = integrator.compute(integrand) + return result.value, result.error + + def _compute_logpdf(self, x: float, integrator: Integrator=LogRQMCIntegrator()) -> tuple[float, float]: + gamma = getattr(self.params, 'gamma', 1) + d = (x - self.params.alpha) ** 2 / gamma ** 2 + + def integrand(u: float) -> float: + return self._log_integrand_func(u, d, gamma) - def compute_cdf(self, x: float, params: dict) -> tuple[float, float]: - parametric_norm = norm(0, self.params.gamma if isinstance(self.params, _NVMClassicDataCollector) else 1) - rqmc = RQMC( - lambda u: parametric_norm.cdf((x - self.params.alpha) / np.sqrt(self.params.distribution.ppf(u))), **params - ) - return rqmc() + result = integrator.compute(integrand) + return result.value, result.error @lru_cache() def _integrand_func(self, u: float, d: float, gamma: float) -> float: ppf = self.params.distribution.ppf(u) - return (1 / np.sqrt(np.pi * 2 * ppf * np.abs(gamma**2))) * np.exp(-1 * d / (2 * ppf)) + return (1 / np.sqrt(2 * np.pi * ppf * np.abs(gamma ** 2))) * np.exp(-d / (2 * ppf)) - def _log_integrand_func(self, u: float, d: float, gamma: float | int | np.int64) -> float: + def _log_integrand_func(self, u: float, d: float, gamma: float) -> float: 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]: - 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() - - def compute_logpdf(self, x: float, params: dict) -> tuple[float, float]: - gamma = self.params.gamma if isinstance(self.params, _NVMClassicDataCollector) else 1 - d = (x - self.params.alpha) ** 2 / gamma**2 - log_rqmc = LogRQMC(lambda u: self._log_integrand_func(u, d, gamma), **params) - return log_rqmc() + return -(ppf * np.log(2 * np.pi * ppf * gamma ** 2) + d) / (2 * ppf) diff --git a/src/procedures/__init__.py b/src/procedures/__init__.py new file mode 100644 index 0000000..00d3e1b --- /dev/null +++ b/src/procedures/__init__.py @@ -0,0 +1,46 @@ +from src.procedures.semiparametric.nm_semiparametric.g_estimation_convolution import ( + NMEstimationDensityInvMT, +) +from src.procedures.semiparametric.nm_semiparametric.sigma_estimation_eigenvalue_based import ( + NMEstimationSigmaEigenvals, +) +from src.procedures.semiparametric.nm_semiparametric.sigma_estimation_empirical import ( + NMEstimationSigmaLaplace, +) +from src.procedures.semiparametric.nv_semiparametric.g_estimation_given_mu import ( + NVEstimationDensityInvMT, +) +from src.procedures.semiparametric.nv_semiparametric.g_estimation_post_widder import ( + NVEstimationDensityPW, +) +from src.procedures.semiparametric.nvm_semiparametric.g_estimation_given_mu import ( + NMVEstimationDensityInvMTquad, +) +from src.procedures.semiparametric.nvm_semiparametric.g_estimation_given_mu_rqmc_based import ( + SemiParametricGEstimationGivenMuRQMCBased, +) +from src.procedures.semiparametric.nvm_semiparametric.g_estimation_post_widder import ( + NMVEstimationDensityPW, +) +from src.procedures.semiparametric.nvm_semiparametric.mu_estimation import (NMVEstimationMu) +from src.register.algorithm_purpose import AlgorithmPurpose +from src.register.register import Registry + +ALGORITHM_REGISTRY: Registry = Registry() +ALGORITHM_REGISTRY.register("mu_estimation", AlgorithmPurpose.NMV_SEMIPARAMETRIC)(NMVEstimationMu) +ALGORITHM_REGISTRY.register("density_estim_inv_mellin_quad_int", AlgorithmPurpose.NMV_SEMIPARAMETRIC)( + NMVEstimationDensityInvMTquad +) +ALGORITHM_REGISTRY.register("density_estim_deconv", AlgorithmPurpose.NM_SEMIPARAMETRIC)(NMEstimationDensityInvMT) +ALGORITHM_REGISTRY.register("sigma_estimation_eigenvalue", AlgorithmPurpose.NM_SEMIPARAMETRIC)( + NMEstimationSigmaEigenvals +) +ALGORITHM_REGISTRY.register("sigma_estimation_laplace", AlgorithmPurpose.NM_SEMIPARAMETRIC)( + NMEstimationSigmaLaplace +) +ALGORITHM_REGISTRY.register("density_estim_inv_mellin_rqmc_int", AlgorithmPurpose.NMV_SEMIPARAMETRIC)( + SemiParametricGEstimationGivenMuRQMCBased +) +ALGORITHM_REGISTRY.register("density_estim_post_widder", AlgorithmPurpose.NMV_SEMIPARAMETRIC)( + NMVEstimationDensityPW +) diff --git a/src/algorithms/param_algorithms/__init__.py b/src/procedures/parametric/__init__.py similarity index 100% rename from src/algorithms/param_algorithms/__init__.py rename to src/procedures/parametric/__init__.py diff --git a/src/algorithms/param_algorithms/nm_param_algorithms/__init__.py b/src/procedures/parametric/nm_parametric/__init__.py similarity index 100% rename from src/algorithms/param_algorithms/nm_param_algorithms/__init__.py rename to src/procedures/parametric/nm_parametric/__init__.py diff --git a/src/algorithms/param_algorithms/nv_param_algorithms/__init__.py b/src/procedures/parametric/nv_parametric/__init__.py similarity index 100% rename from src/algorithms/param_algorithms/nv_param_algorithms/__init__.py rename to src/procedures/parametric/nv_parametric/__init__.py diff --git a/src/algorithms/param_algorithms/nvm_param_algorithms/__init__.py b/src/procedures/parametric/nvm_parametric/__init__.py similarity index 100% rename from src/algorithms/param_algorithms/nvm_param_algorithms/__init__.py rename to src/procedures/parametric/nvm_parametric/__init__.py diff --git a/src/algorithms/semiparam_algorithms/__init__.py b/src/procedures/semiparametric/__init__.py similarity index 100% rename from src/algorithms/semiparam_algorithms/__init__.py rename to src/procedures/semiparametric/__init__.py diff --git a/src/algorithms/semiparam_algorithms/nm_semi_param_algorithms/__init__.py b/src/procedures/semiparametric/nm_semiparametric/__init__.py similarity index 100% rename from src/algorithms/semiparam_algorithms/nm_semi_param_algorithms/__init__.py rename to src/procedures/semiparametric/nm_semiparametric/__init__.py diff --git a/src/algorithms/semiparam_algorithms/nm_semi_param_algorithms/g_estimation_convolution.py b/src/procedures/semiparametric/nm_semiparametric/g_estimation_convolution.py similarity index 91% rename from src/algorithms/semiparam_algorithms/nm_semi_param_algorithms/g_estimation_convolution.py rename to src/procedures/semiparametric/nm_semiparametric/g_estimation_convolution.py index 3feb6cf..83de09f 100644 --- a/src/algorithms/semiparam_algorithms/nm_semi_param_algorithms/g_estimation_convolution.py +++ b/src/procedures/semiparametric/nm_semiparametric/g_estimation_convolution.py @@ -11,8 +11,7 @@ BOHMAN_DELTA_DEFAULT_VALUE: float = 0.0001 X_DATA_DEFAULT_VALUE: List[float] = [1.0] - -class NMSemiParametricGEstimation: +class NMEstimationDensityInvMT: """Estimation of mixing density function g (xi density function) of NM mixture represented in canonical form Y = xi + sigma*N. @@ -34,8 +33,8 @@ class ParamsAnnotation(TypedDict, total=False): bohman_n: int bohman_delta: float - def __init__(self, sample: Optional[_typing.NDArray[np.float64]] = None, **kwargs: Unpack[ParamsAnnotation]): - self.sample: _typing.NDArray[np.float64] = np.array([]) if sample is None else sample + def __init__(self, sample: Optional[np.ndarray] = None, **kwargs: Unpack[ParamsAnnotation]): + self.sample: np.ndarray = np.array([]) if sample is None else sample self.n: int = len(self.sample) ( self.x_data, @@ -93,7 +92,7 @@ def cdf(self, X: np.ndarray) -> np.ndarray: return F_x.real - def algorithm(self, sample: _typing.NDArray[np.float64]) -> EstimateResult: + def compute(self, sample: np.ndarray) -> EstimateResult: inv = self.BohmanA(N=self.bohman_n, delta=self.bohman_delta) inv.fit(self.characteristic_function_xi) x_data_array = np.array(self.x_data, dtype=np.float64) diff --git a/src/algorithms/semiparam_algorithms/nm_semi_param_algorithms/sigma_estimation_eigenvalue_based.py b/src/procedures/semiparametric/nm_semiparametric/sigma_estimation_eigenvalue_based.py similarity index 97% rename from src/algorithms/semiparam_algorithms/nm_semi_param_algorithms/sigma_estimation_eigenvalue_based.py rename to src/procedures/semiparametric/nm_semiparametric/sigma_estimation_eigenvalue_based.py index 1148610..b9da043 100644 --- a/src/algorithms/semiparam_algorithms/nm_semi_param_algorithms/sigma_estimation_eigenvalue_based.py +++ b/src/procedures/semiparametric/nm_semiparametric/sigma_estimation_eigenvalue_based.py @@ -14,7 +14,7 @@ PARAMETER_KEYS = ["l", "k", "eps", "search_area", "search_density"] -class SemiParametricMeanSigmaEstimationEigenvalueBased: +class NMEstimationSigmaEigenvals: """Estimation of sigma parameter of NM mixture represented in canonical form Y = xi + sigma*N. Args: @@ -86,7 +86,7 @@ def _build_matrix(self, tau: float) -> np.ndarray: matrix[i, j] = self._alpha(t[i] - t[j], tau) return matrix - def algorithm(self, sample: np.ndarray) -> EstimateResult: + def compute(self, sample: np.ndarray) -> EstimateResult: """Estimate sigma. Args: diff --git a/src/algorithms/semiparam_algorithms/nm_semi_param_algorithms/sigma_estimation_empirical.py b/src/procedures/semiparametric/nm_semiparametric/sigma_estimation_empirical.py similarity index 94% rename from src/algorithms/semiparam_algorithms/nm_semi_param_algorithms/sigma_estimation_empirical.py rename to src/procedures/semiparametric/nm_semiparametric/sigma_estimation_empirical.py index 0ce471c..22a69af 100644 --- a/src/algorithms/semiparam_algorithms/nm_semi_param_algorithms/sigma_estimation_empirical.py +++ b/src/procedures/semiparametric/nm_semiparametric/sigma_estimation_empirical.py @@ -9,7 +9,7 @@ PARAMETER_KEYS = ["t"] -class SemiParametricMeanSigmaEstimationEmpirical: +class NMEstimationSigmaLaplace: """Estimation of sigma parameter of NM mixture represented in canonical form Y = xi + sigma*N. Args: @@ -51,7 +51,7 @@ def _validate_kwargs(**kwargs: ParamsAnnotation) -> float: raise ValueError("Expected a positive float as parameter `t`.") return t - def algorithm(self, sample: np.ndarray) -> EstimateResult: + def compute(self, sample: np.ndarray) -> EstimateResult: """Estimate the sigma parameter. Args: diff --git a/src/algorithms/semiparam_algorithms/nv_semi_param_algorithms/__init__.py b/src/procedures/semiparametric/nv_semiparametric/__init__.py similarity index 100% rename from src/algorithms/semiparam_algorithms/nv_semi_param_algorithms/__init__.py rename to src/procedures/semiparametric/nv_semiparametric/__init__.py diff --git a/src/algorithms/semiparam_algorithms/nv_semi_param_algorithms/g_estimation_given_mu.py b/src/procedures/semiparametric/nv_semiparametric/g_estimation_given_mu.py similarity index 95% rename from src/algorithms/semiparam_algorithms/nv_semi_param_algorithms/g_estimation_given_mu.py rename to src/procedures/semiparametric/nv_semiparametric/g_estimation_given_mu.py index 7c810cc..768d54c 100644 --- a/src/algorithms/semiparam_algorithms/nv_semi_param_algorithms/g_estimation_given_mu.py +++ b/src/procedures/semiparametric/nv_semiparametric/g_estimation_given_mu.py @@ -30,7 +30,7 @@ def v_sequence_default_value(n: float) -> float: INTEGRATION_LIMIT_DEFAULT_VALUE: int = 50 -class SemiParametricNVEstimation: +class NVEstimationDensityInvMT: """Estimation of mixing density function g (xi density function) of NV mixture represented in canonical form Y = alpha + sqrt(xi)*N, where alpha = 0 and mu = 0. @@ -49,13 +49,13 @@ class ParamsAnnotation(TypedDict, total=False): integration_tolerance: float integration_limit: int - def __init__(self, sample: Optional[_typing.NDArray[np.float64]] = None, **kwargs: Unpack[ParamsAnnotation]): + def __init__(self, sample: Optional[np.ndarray] = None, **kwargs: Unpack[ParamsAnnotation]): self.x_powers: Dict[float, np.ndarray] = {} self.second_u_integrals: np.ndarray self.first_u_integrals: np.ndarray self.gamma_grid: np.ndarray self.v_grid: np.ndarray - self.sample: _typing.NDArray[np.float64] = np.array([]) if sample is None else sample + self.sample: np.ndarray = np.array([]) if sample is None else sample self.n: int = len(self.sample) ( self.gmm, @@ -173,6 +173,6 @@ def compute_integrals_for_x(self, x: float) -> float: total = np.sum(first_integral + second_integral) / self.denominator return max(0.0, total.real) - def algorithm(self, sample: _typing.NDArray[np.float64]) -> EstimateResult: + def compute(self, sample: np.ndarray) -> EstimateResult: y_data = [self.compute_integrals_for_x(x) for x in self.x_data] return EstimateResult(list_value=y_data, success=True) diff --git a/src/algorithms/semiparam_algorithms/nv_semi_param_algorithms/g_estimation_post_widder.py b/src/procedures/semiparametric/nv_semiparametric/g_estimation_post_widder.py similarity index 96% rename from src/algorithms/semiparam_algorithms/nv_semi_param_algorithms/g_estimation_post_widder.py rename to src/procedures/semiparametric/nv_semiparametric/g_estimation_post_widder.py index 0c681d3..14c6deb 100644 --- a/src/algorithms/semiparam_algorithms/nv_semi_param_algorithms/g_estimation_post_widder.py +++ b/src/procedures/semiparametric/nv_semiparametric/g_estimation_post_widder.py @@ -13,7 +13,7 @@ X_DATA_DEFAULT_VALUE = [1.0] -class NVSemiParametricGEstimationPostWidder: +class NVEstimationDensityPW: """Estimation of mixing density function g (xi density function) of NVM mixture represented in classical form Y = alpha + sigma*sqrt(xi)*N, where alpha = 0 and mu, sigma are given. @@ -77,7 +77,7 @@ def p_x_estimation(self, x: float) -> float: ) return result.real - def algorithm(self, sample: np._typing.NDArray) -> EstimateResult: + def compute(self, sample: np.ndarray) -> EstimateResult: """Estimate g(x) Args: diff --git a/src/algorithms/semiparam_algorithms/nvm_semi_param_algorithms/__init__.py b/src/procedures/semiparametric/nvm_semiparametric/__init__.py similarity index 100% rename from src/algorithms/semiparam_algorithms/nvm_semi_param_algorithms/__init__.py rename to src/procedures/semiparametric/nvm_semiparametric/__init__.py diff --git a/src/algorithms/semiparam_algorithms/nvm_semi_param_algorithms/g_estimation_given_mu.py b/src/procedures/semiparametric/nvm_semiparametric/g_estimation_given_mu.py similarity index 95% rename from src/algorithms/semiparam_algorithms/nvm_semi_param_algorithms/g_estimation_given_mu.py rename to src/procedures/semiparametric/nvm_semiparametric/g_estimation_given_mu.py index cbda86e..b84363c 100644 --- a/src/algorithms/semiparam_algorithms/nvm_semi_param_algorithms/g_estimation_given_mu.py +++ b/src/procedures/semiparametric/nvm_semiparametric/g_estimation_given_mu.py @@ -31,7 +31,7 @@ def v_sequence_default_value(n: float) -> float: INTEGRATION_LIMIT_DEFAULT_VALUE: int = 50 -class SemiParametricGEstimationGivenMu: +class NMVEstimationDensityInvMTquad: """Estimation of mixing density function g (xi density function) of NVM mixture represented in canonical form Y = alpha + mu*xi + sqrt(xi)*N, where alpha = 0 and mu is given. @@ -51,13 +51,13 @@ class ParamsAnnotation(TypedDict, total=False): integration_tolerance: float integration_limit: int - def __init__(self, sample: Optional[_typing.NDArray[np.float64]] = None, **kwargs: Unpack[ParamsAnnotation]): + def __init__(self, sample: Optional[np.ndarray] = None, **kwargs: Unpack[ParamsAnnotation]): self.x_powers: Dict[float, np.ndarray] = {} self.second_u_integrals: np.ndarray self.first_u_integrals: np.ndarray self.gamma_grid: np.ndarray self.v_grid: np.ndarray - self.sample: _typing.NDArray[np.float64] = np.array([]) if sample is None else sample + self.sample: np.ndarray= np.array([]) if sample is None else sample self.n: int = len(self.sample) ( self.mu, @@ -179,6 +179,6 @@ def compute_integrals_for_x(self, x: float) -> float: total = np.sum(first_integral + second_integral) / self.denominator return max(0.0, total.real) - def algorithm(self, sample: _typing.NDArray[np.float64]) -> EstimateResult: + def compute(self, sample: np.ndarray) -> EstimateResult: y_data = [self.compute_integrals_for_x(x) for x in self.x_data] return EstimateResult(list_value=y_data, success=True) diff --git a/src/algorithms/semiparam_algorithms/nvm_semi_param_algorithms/g_estimation_given_mu_rqmc_based.py b/src/procedures/semiparametric/nvm_semiparametric/g_estimation_given_mu_rqmc_based.py similarity index 90% rename from src/algorithms/semiparam_algorithms/nvm_semi_param_algorithms/g_estimation_given_mu_rqmc_based.py rename to src/procedures/semiparametric/nvm_semiparametric/g_estimation_given_mu_rqmc_based.py index 2f3da81..2bf1687 100644 --- a/src/algorithms/semiparam_algorithms/nvm_semi_param_algorithms/g_estimation_given_mu_rqmc_based.py +++ b/src/procedures/semiparametric/nvm_semiparametric/g_estimation_given_mu_rqmc_based.py @@ -3,11 +3,11 @@ from typing import Callable, Dict, List, Optional, TypedDict, Unpack import numpy as np -from numpy import _typing from scipy.integrate import quad_vec from scipy.special import gamma -from src.algorithms.support_algorithms.rqmc import RQMC +from src.procedures.support.integrator import Integrator +from src.procedures.support.rqmc import RQMCIntegrator from src.estimators.estimate_result import EstimateResult MU_DEFAULT_VALUE = 1.0 @@ -52,13 +52,13 @@ class ParamsAnnotation(TypedDict, total=False): integration_tolerance: float integration_limit: int - def __init__(self, sample: Optional[_typing.NDArray[np.float64]] = None, **kwargs: Unpack[ParamsAnnotation]): + def __init__(self, sample: Optional[np.ndarray] = None, **kwargs: Unpack[ParamsAnnotation]): self.x_powers: Dict[float, np.ndarray] = {} self.second_u_integrals: np.ndarray self.first_u_integrals: np.ndarray self.gamma_grid: np.ndarray self.v_grid: np.ndarray - self.sample: _typing.NDArray[np.float64] = np.array([]) if sample is None else sample + self.sample: np.ndarray = np.array([]) if sample is None else sample self.n: int = len(self.sample) ( self.mu, @@ -68,7 +68,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.integration_limit ) = self._validate_kwargs(self.n, **kwargs) self.denominator: float = 2 * math.pi * self.n self.precompute_gamma_grid() @@ -164,13 +164,14 @@ def second_v_integrand(self, v: float, x: float) -> np.ndarray: def compute_integrals_for_x(self, x: float) -> 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 = RQMCIntegrator() + first_integral = integrator.compute(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(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: + def algorithm(self, sample: np.ndarray) -> EstimateResult: y_data = [self.compute_integrals_for_x(x) for x in self.x_data] return EstimateResult(list_value=y_data, success=True) diff --git a/src/algorithms/semiparam_algorithms/nvm_semi_param_algorithms/g_estimation_post_widder.py b/src/procedures/semiparametric/nvm_semiparametric/g_estimation_post_widder.py similarity index 96% rename from src/algorithms/semiparam_algorithms/nvm_semi_param_algorithms/g_estimation_post_widder.py rename to src/procedures/semiparametric/nvm_semiparametric/g_estimation_post_widder.py index 50c518e..2d219b0 100644 --- a/src/algorithms/semiparam_algorithms/nvm_semi_param_algorithms/g_estimation_post_widder.py +++ b/src/procedures/semiparametric/nvm_semiparametric/g_estimation_post_widder.py @@ -14,7 +14,7 @@ X_DATA_DEFAULT_VALUE = [1.0] -class SemiParametricGEstimationPostWidder: +class NMVEstimationDensityPW: """Estimation of mixing density function g (xi density function) of NVM mixture represented in classical form Y = alpha + mu*xi + sigma*sqrt(xi)*N, where alpha = 0 and mu, sigma are given. @@ -80,7 +80,7 @@ def p_x_estimation(self, x: float) -> float: ) return result.real - def algorithm(self, sample: np._typing.NDArray) -> EstimateResult: + def compute(self, sample: np._typing.NDArray) -> EstimateResult: """Estimate g(x) Args: diff --git a/src/algorithms/semiparam_algorithms/nvm_semi_param_algorithms/mu_estimation.py b/src/procedures/semiparametric/nvm_semiparametric/mu_estimation.py similarity index 96% rename from src/algorithms/semiparam_algorithms/nvm_semi_param_algorithms/mu_estimation.py rename to src/procedures/semiparametric/nvm_semiparametric/mu_estimation.py index 9e2dad0..3ca6292 100644 --- a/src/algorithms/semiparam_algorithms/nvm_semi_param_algorithms/mu_estimation.py +++ b/src/procedures/semiparametric/nvm_semiparametric/mu_estimation.py @@ -14,7 +14,7 @@ PARAMETER_KEYS = ["m", "tolerance", "omega", "max_iterations"] -class SemiParametricMuEstimation: +class NMVEstimationMu: """Estimation of mu parameter of NVM mixture represented in canonical form Y = alpha + mu*xi + sqrt(xi)*N, where alpha = 0 @@ -84,7 +84,7 @@ def __w(self, p: float, sample: np._typing.NDArray) -> float: y += e * self.omega(x) return y - def algorithm(self, sample: np._typing.NDArray) -> EstimateResult: + def compute(self, sample: np.ndarray) -> EstimateResult: """Root of this function is an estimation of mu Args: @@ -96,7 +96,7 @@ def algorithm(self, sample: np._typing.NDArray) -> EstimateResult: if self.__w(0, sample) == 0: return EstimateResult(value=0, success=True) if self.__w(0, sample) > 0: - second_result = self.algorithm(-1 * sample) + second_result = self.compute(-1 * sample) return EstimateResult(value=-1 * second_result.value, success=second_result.success) if self.__w(self.m, sample) < 0: return EstimateResult(value=self.m, success=False) diff --git a/src/procedures/semiparametric/statistical_procedure.py b/src/procedures/semiparametric/statistical_procedure.py new file mode 100644 index 0000000..228c82c --- /dev/null +++ b/src/procedures/semiparametric/statistical_procedure.py @@ -0,0 +1,12 @@ +from typing import Protocol +import numpy as np + +from src.estimators.estimate_result import EstimateResult + +class StatisticalProcedure(Protocol): + """Protocol for statistical procedures that compute estimates from data. + + """ + def compute(self, sample: np.ndarray) -> EstimateResult: + pass + diff --git a/src/procedures/support/integrator.py b/src/procedures/support/integrator.py new file mode 100644 index 0000000..8ee0fb7 --- /dev/null +++ b/src/procedures/support/integrator.py @@ -0,0 +1,20 @@ +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 __init__(self) -> None: + ... + + def compute(self, func: Callable) -> IntegrationResult: + ... diff --git a/src/procedures/support/log_rqmc.py b/src/procedures/support/log_rqmc.py new file mode 100644 index 0000000..af267b0 --- /dev/null +++ b/src/procedures/support/log_rqmc.py @@ -0,0 +1,81 @@ +from typing import Callable +import numpy as np +import numpy._typing as tpg +import scipy +from src.procedures.support.integrator import IntegrationResult +from src.procedures.support.rqmc import RQMC + + +class LogRQMCIntegrator: + """Log-sum-exp stabilized randomized Quasi-Monte Carlo integrator.""" + + 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 + + @staticmethod + def lse(args: tpg.NDArray) -> float: + """ + Log-Sum-Exp stabilization helper. + """ + max_value = max(args) + return max_value + np.log(np.sum(np.exp(args - max_value))) + + def _independent_estimator(self, values: tpg.NDArray) -> float: + vfunc = np.vectorize(self.func) + return -np.log(len(values)) + self.lse(vfunc(values)) + + def _estimator(self, random_matrix: tpg.NDArray) -> tuple[float, tpg.NDArray]: + values = np.array(list(map(self._independent_estimator, random_matrix))) + return -np.log(self.count) + self.lse(values), values + + def _update_independent_estimator( + self, i: int, old_value: float, new_values: tpg.NDArray + ) -> float: + return -np.log(i + 1) + self.lse( + np.array(i * [old_value] + [self._independent_estimator(new_values[i * self.base_n :])]) + ) + + def _update( + self, j: int, old_values: tpg.NDArray, random_matrix: tpg.NDArray + ) -> tuple[float, tpg.NDArray]: + values = [] + for idx in range(self.count): + old_val, new_vals = old_values[idx], random_matrix[idx] + values.append(self._update_independent_estimator(j, old_val, new_vals)) + np_values = np.array(values) + return -np.log(self.count) + self.lse(np_values), np_values + + def compute(self, func: Callable) -> IntegrationResult: + """Compute integral of `func` over [0,1] with log-RQMC.""" + # Assign the integrand + self.func = func + # Instantiate underlying RQMC with same parameters + rqmc_inst = RQMC( + func, + error_tolerance=self.error_tolerance, + count=self.count, + base_n=self.base_n, + i_max=self.i_max, + a=self.a, + ) + # __call__ returns (approximation, error) + approximation, error = rqmc_inst() + return IntegrationResult(approximation, error, None) + + +if __name__ == "__main__": + # Sample usage + integrator = LogRQMCIntegrator(i_max=100) + result = integrator.compute(lambda x: x**3 - x**2 + 1) + print(result) diff --git a/src/procedures/support/quad_integrator.py b/src/procedures/support/quad_integrator.py new file mode 100644 index 0000000..8401887 --- /dev/null +++ b/src/procedures/support/quad_integrator.py @@ -0,0 +1,60 @@ +from typing import Callable, Any, Sequence + +from scipy.integrate import quad +from src.procedures.support.integrator import IntegrationResult + +class QuadIntegrator: + + 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 + + Returns: moment approximation and error tolerance + """ + + verbose = self.params.pop('full_output', False) + result = quad(func, **self.params) + if verbose: + value, error, message = result + else: + value, error = result + message = None + return IntegrationResult(value, error, message) diff --git a/src/algorithms/support_algorithms/rqmc.py b/src/procedures/support/rqmc.py similarity index 83% rename from src/algorithms/support_algorithms/rqmc.py rename to src/procedures/support/rqmc.py index 4b10c9d..4254615 100644 --- a/src/algorithms/support_algorithms/rqmc.py +++ b/src/procedures/support/rqmc.py @@ -5,6 +5,8 @@ import scipy from numba import njit +from src.procedures.support.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,51 @@ def __call__(self) -> tuple[float, float]: """ return self.rqmc() + + +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(self, func: Callable) -> IntegrationResult: + """ + Compute integral via RQMC integrator + + Args: + func: integrated function + + 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]) + diff --git a/tests/generators/nm_generator/test_mixing_normal.py b/tests/generators/nm_generator/test_mixing_normal.py index ea85692..49af0ae 100644 --- a/tests/generators/nm_generator/test_mixing_normal.py +++ b/tests/generators/nm_generator/test_mixing_normal.py @@ -16,7 +16,7 @@ class TestMixingNormal: ) def test_classic_generate_variance_0(self, mixing_variance: float, expected_variance: float) -> None: mixture = NormalMeanMixtures("classical", alpha=0, beta=mixing_variance**0.5, gamma=1, distribution=norm) - sample = self.generator.classical_generate(mixture, self.test_mixture_size) + sample = self.generator.generate(mixture, self.test_mixture_size) actual_variance = ndimage.variance(sample) assert actual_variance == pytest.approx(expected_variance, 0.1) @@ -24,7 +24,7 @@ def test_classic_generate_variance_0(self, mixing_variance: float, expected_vari def test_classic_generate_variance_1(self, beta: float) -> None: expected_variance = beta**2 + 1 mixture = NormalMeanMixtures("classical", alpha=0, beta=beta, gamma=1, distribution=norm) - sample = self.generator.classical_generate(mixture, self.test_mixture_size) + sample = self.generator.generate(mixture, self.test_mixture_size) actual_variance = ndimage.variance(sample) assert actual_variance == pytest.approx(expected_variance, 0.1) @@ -32,7 +32,7 @@ def test_classic_generate_variance_1(self, beta: float) -> None: def test_classic_generate_variance_2(self, beta: float, gamma: float) -> None: expected_variance = beta**2 + gamma**2 mixture = NormalMeanMixtures("classical", alpha=0, beta=beta, gamma=gamma, distribution=norm) - sample = self.generator.classical_generate(mixture, self.test_mixture_size) + sample = self.generator.generate(mixture, self.test_mixture_size) actual_variance = ndimage.variance(sample) assert actual_variance == pytest.approx(expected_variance, 0.1) @@ -40,14 +40,14 @@ def test_classic_generate_variance_2(self, beta: float, gamma: float) -> None: def test_classic_generate_mean(self, beta: float, gamma: float) -> None: expected_mean = 0 mixture = NormalMeanMixtures("classical", alpha=0, beta=beta, gamma=gamma, distribution=norm) - sample = self.generator.classical_generate(mixture, self.test_mixture_size) + sample = self.generator.generate(mixture, self.test_mixture_size) actual_mean = np.mean(np.array(sample)) assert abs(actual_mean - expected_mean) < 1 @pytest.mark.parametrize("expected_size", np.random.randint(0, 100, size=50)) def test_classic_generate_size(self, expected_size: int) -> None: mixture = NormalMeanMixtures("classical", alpha=0, beta=1, gamma=1, distribution=norm) - sample = self.generator.classical_generate(mixture, expected_size) + sample = self.generator.generate(mixture, expected_size) actual_size = np.size(sample) assert actual_size == expected_size @@ -56,7 +56,7 @@ def test_classic_generate_size(self, expected_size: int) -> None: ) def test_canonical_generate_variance_0(self, mixing_variance: float, expected_variance: float) -> None: mixture = NormalMeanMixtures("canonical", sigma=1, distribution=norm(0, mixing_variance**0.5)) - sample = self.generator.canonical_generate(mixture, self.test_mixture_size) + sample = self.generator.generate(mixture, self.test_mixture_size) actual_variance = ndimage.variance(sample) assert actual_variance == pytest.approx(expected_variance, 0.1) @@ -64,7 +64,7 @@ def test_canonical_generate_variance_0(self, mixing_variance: float, expected_va def test_canonical_generate_variance_1(self, sigma: float) -> None: expected_variance = sigma**2 + 1 mixture = NormalMeanMixtures("canonical", sigma=sigma, distribution=norm) - sample = self.generator.canonical_generate(mixture, self.test_mixture_size) + sample = self.generator.generate(mixture, self.test_mixture_size) actual_variance = ndimage.variance(sample) assert actual_variance == pytest.approx(expected_variance, 0.1) @@ -72,7 +72,7 @@ def test_canonical_generate_variance_1(self, sigma: float) -> None: def test_canonical_generate_variance_2(self, mixing_variance: float, sigma: float) -> None: expected_variance = mixing_variance + sigma**2 mixture = NormalMeanMixtures("canonical", sigma=sigma, distribution=norm(0, mixing_variance**0.5)) - sample = self.generator.canonical_generate(mixture, self.test_mixture_size) + sample = self.generator.generate(mixture, self.test_mixture_size) actual_variance = ndimage.variance(sample) assert actual_variance == pytest.approx(expected_variance, 0.1) @@ -80,13 +80,13 @@ def test_canonical_generate_variance_2(self, mixing_variance: float, sigma: floa def test_canonical_generate_mean(self, sigma: float) -> None: expected_mean = 0 mixture = NormalMeanMixtures("canonical", sigma=sigma, distribution=norm) - sample = self.generator.canonical_generate(mixture, self.test_mixture_size) + sample = self.generator.generate(mixture, self.test_mixture_size) actual_mean = np.mean(np.array(sample)) assert abs(actual_mean - expected_mean) < 1 @pytest.mark.parametrize("expected_size", [*np.random.randint(0, 100, size=50), 0, 1, 1000000]) def test_canonical_generate_size(self, expected_size: int) -> None: mixture = NormalMeanMixtures("canonical", sigma=1, distribution=norm) - sample = self.generator.canonical_generate(mixture, expected_size) + sample = self.generator.generate(mixture, expected_size) actual_size = np.size(sample) assert actual_size == expected_size diff --git a/tests/algorithms/__init__.py b/tests/procedures/__init__.py similarity index 100% rename from tests/algorithms/__init__.py rename to tests/procedures/__init__.py diff --git a/tests/algorithms/nm_algorithms/semiparametric_sigma_estimation/__init__.py b/tests/procedures/nm_procedures/semiparametric_sigma_estimation/__init__.py similarity index 100% rename from tests/algorithms/nm_algorithms/semiparametric_sigma_estimation/__init__.py rename to tests/procedures/nm_procedures/semiparametric_sigma_estimation/__init__.py diff --git a/tests/algorithms/nm_algorithms/semiparametric_sigma_estimation/test_sigma_estimation_eigenvalue_based.py b/tests/procedures/nm_procedures/semiparametric_sigma_estimation/test_sigma_estimation_eigenvalue_based.py similarity index 86% rename from tests/algorithms/nm_algorithms/semiparametric_sigma_estimation/test_sigma_estimation_eigenvalue_based.py rename to tests/procedures/nm_procedures/semiparametric_sigma_estimation/test_sigma_estimation_eigenvalue_based.py index eb00242..e3b7d64 100644 --- a/tests/algorithms/nm_algorithms/semiparametric_sigma_estimation/test_sigma_estimation_eigenvalue_based.py +++ b/tests/procedures/nm_procedures/semiparametric_sigma_estimation/test_sigma_estimation_eigenvalue_based.py @@ -4,7 +4,7 @@ import pytest from scipy.stats import expon, uniform -from src.estimators.semiparametric.nm_semiparametric_estimator import NMSemiParametricEstimator +from src.estimators.semiparametric.nm_semiparametric_estimator import NMSemiparEstim from src.generators.nm_generator import NMGenerator from src.mixtures.nm_mixture import NormalMeanMixtures @@ -28,8 +28,8 @@ def test_sigma_estimation_eigenvalue_based_expon_single( * (np.log(np.log(sample_len))) ** 0.6 ) - estimator = NMSemiParametricEstimator( - "sigma_estimation_eigenvalue_based", {"k": k, "l": l, "eps": eps, "search_area": search_area} + estimator = NMSemiparEstim( + "sigma_estimation_eigenvalue", {"k": k, "l": l, "eps": eps, "search_area": search_area} ) est = estimator.estimate(sample) print(est.value) @@ -54,8 +54,8 @@ def test_sigma_estimation_eigenvalue_based_expon_1( * (np.log(np.log(sample_len))) ** 0.6 ) - estimator = NMSemiParametricEstimator( - "sigma_estimation_eigenvalue_based", {"k": k, "l": l, "eps": eps, "search_area": search_area} + estimator = NMSemiparEstim( + "sigma_estimation_eigenvalue", {"k": k, "l": l, "eps": eps, "search_area": search_area} ) est = estimator.estimate(sample) print(est.value) @@ -80,8 +80,8 @@ def test_sigma_estimation_eigenvalue_based_expon_2( * (np.log(np.log(sample_len))) ** 0.6 ) - estimator = NMSemiParametricEstimator( - "sigma_estimation_eigenvalue_based", {"k": k, "l": l, "eps": eps, "search_area": search_area} + estimator = NMSemiparEstim( + "sigma_estimation_eigenvalue", {"k": k, "l": l, "eps": eps, "search_area": search_area} ) est = estimator.estimate(sample) print(est.value) @@ -106,8 +106,8 @@ def test_sigma_estimation_eigenvalue_based_expon_3( * (np.log(np.log(sample_len))) ** 0.6 ) - estimator = NMSemiParametricEstimator( - "sigma_estimation_eigenvalue_based", {"k": k, "l": l, "eps": eps, "search_area": search_area} + estimator = NMSemiparEstim( + "sigma_estimation_eigenvalue", {"k": k, "l": l, "eps": eps, "search_area": search_area} ) est = estimator.estimate(sample) print(est.value) @@ -131,8 +131,8 @@ def test_sigma_estimation_eigenvalue_based_uniform( * (np.log(np.log(sample_len))) ** 0.6 ) - estimator = NMSemiParametricEstimator( - "sigma_estimation_eigenvalue_based", {"k": k, "l": l, "eps": eps, "search_area": search_area} + estimator = NMSemiparEstim( + "sigma_estimation_eigenvalue", {"k": k, "l": l, "eps": eps, "search_area": search_area} ) est = estimator.estimate(sample) print(est.value) diff --git a/tests/algorithms/nm_algorithms/semiparametric_sigma_estimation/test_sigma_estimation_empirical.py b/tests/procedures/nm_procedures/semiparametric_sigma_estimation/test_sigma_estimation_empirical.py similarity index 91% rename from tests/algorithms/nm_algorithms/semiparametric_sigma_estimation/test_sigma_estimation_empirical.py rename to tests/procedures/nm_procedures/semiparametric_sigma_estimation/test_sigma_estimation_empirical.py index daf5a7f..550978f 100644 --- a/tests/algorithms/nm_algorithms/semiparametric_sigma_estimation/test_sigma_estimation_empirical.py +++ b/tests/procedures/nm_procedures/semiparametric_sigma_estimation/test_sigma_estimation_empirical.py @@ -3,7 +3,7 @@ import pytest from scipy.stats import expon, uniform -from src.estimators.semiparametric.nm_semiparametric_estimator import NMSemiParametricEstimator +from src.estimators.semiparametric.nm_semiparametric_estimator import NMSemiparEstim from src.generators.nm_generator import NMGenerator from src.mixtures.nm_mixture import NormalMeanMixtures @@ -20,7 +20,7 @@ def test_sigma_estimation_empirical_expon(self, real_sigma: float, sample_len: i answer_list = [] for alpha in [round(x, 5) for x in [i * 0.0001 for i in range(1, 10000)]]: t = math.sqrt(alpha * math.log(sample_len)) / (2 * search_area) - estimator = NMSemiParametricEstimator("sigma_estimation_empirical", {"t": t}) + estimator = NMSemiparEstim("sigma_estimation_laplace", {"t": t}) est = estimator.estimate(sample) left = (est.value**2 - real_sigma**2) ** 0.5 right = ( @@ -39,7 +39,7 @@ def test_sigma_estimation_empirical_uniform(self, real_sigma: float, sample_len: M = 10 for alpha in [round(x, 4) for x in [i * 0.0001 for i in range(1, 10000)]]: t = math.sqrt(alpha * math.log(sample_len)) / (2 * search_area) - estimator = NMSemiParametricEstimator("sigma_estimation_empirical", {"t": t}) + estimator = NMSemiparEstim("sigma_estimation_laplace", {"t": t}) est = estimator.estimate(sample) left = abs(est.value**2 - real_sigma**2) ** 0.5 right = 4 * M * search_area / math.sqrt(alpha * math.log(sample_len)) diff --git a/tests/algorithms/nmv_algorithms/semiparametric_g_estimation/__init__.py b/tests/procedures/nmv_procedures/semiparametric_g_estimation/__init__.py similarity index 100% rename from tests/algorithms/nmv_algorithms/semiparametric_g_estimation/__init__.py rename to tests/procedures/nmv_procedures/semiparametric_g_estimation/__init__.py diff --git a/tests/algorithms/nmv_algorithms/semiparametric_g_estimation/test_g_estimation_given_mu.py b/tests/procedures/nmv_procedures/semiparametric_g_estimation/test_g_estimation_given_mu.py similarity index 84% rename from tests/algorithms/nmv_algorithms/semiparametric_g_estimation/test_g_estimation_given_mu.py rename to tests/procedures/nmv_procedures/semiparametric_g_estimation/test_g_estimation_given_mu.py index d1f40e0..0cc968f 100644 --- a/tests/algorithms/nmv_algorithms/semiparametric_g_estimation/test_g_estimation_given_mu.py +++ b/tests/procedures/nmv_procedures/semiparametric_g_estimation/test_g_estimation_given_mu.py @@ -3,7 +3,7 @@ import pytest from scipy.stats import expon, gamma -from src.estimators.semiparametric.nmv_semiparametric_estimator import NMVSemiParametricEstimator +from src.estimators.semiparametric.nmv_semiparametric_estimator import NMVSemiparEstim from src.generators.nmv_generator import NMVGenerator from src.mixtures.nmv_mixture import NormalMeanVarianceMixtures @@ -18,8 +18,8 @@ def test_g_estimation_expon(self, x_data) -> None: mixture = NormalMeanVarianceMixtures("canonical", alpha=0, mu=given_mu, distribution=expon) sample = NMVGenerator().canonical_generate(mixture, n) - estimator = NMVSemiParametricEstimator( - "g_estimation_given_mu", {"x_data": x_data, "u_value": 7.6, "v_value": 0.9} + estimator = NMVSemiparEstim( + "density_estim_inv_mellin_quad_int", {"x_data": x_data, "u_value": 7.6, "v_value": 0.9} ) est = estimator.estimate(sample) error = 0.0 diff --git a/tests/algorithms/nmv_algorithms/semiparametric_g_estimation/test_post_widder.py b/tests/procedures/nmv_procedures/semiparametric_g_estimation/test_post_widder.py similarity index 85% rename from tests/algorithms/nmv_algorithms/semiparametric_g_estimation/test_post_widder.py rename to tests/procedures/nmv_procedures/semiparametric_g_estimation/test_post_widder.py index 2e8d4b3..e3b268c 100644 --- a/tests/algorithms/nmv_algorithms/semiparametric_g_estimation/test_post_widder.py +++ b/tests/procedures/nmv_procedures/semiparametric_g_estimation/test_post_widder.py @@ -3,7 +3,7 @@ from mpmath import ln from scipy.stats import expon, gamma -from src.estimators.semiparametric.nmv_semiparametric_estimator import NMVSemiParametricEstimator +from src.estimators.semiparametric.nmv_semiparametric_estimator import NMVSemiparEstim from src.generators.nmv_generator import NMVGenerator from src.mixtures.nmv_mixture import NormalMeanVarianceMixtures @@ -20,8 +20,8 @@ def test_post_widder_expon(self, mu, sigma, degree, sample_size) -> None: sample = NMVGenerator().classical_generate(mixture, sample_size) x_data = np.linspace(0.5, 10.0, 30) - estimator = NMVSemiParametricEstimator( - "g_estimation_post_widder", {"x_data": x_data, "mu": mu, "sigma": sigma, "n": degree} + estimator = NMVSemiparEstim( + "density_estim_post_widder", {"x_data": x_data, "mu": mu, "sigma": sigma, "n": degree} ) est = estimator.estimate(sample) est_data = est.list_value @@ -43,8 +43,8 @@ def test_post_widder_gamma(self, mu, sigma, degree, sample_size, a) -> None: sample = NMVGenerator().classical_generate(mixture, sample_size) x_data = np.linspace(0.5, 10.0, 30) - estimator = NMVSemiParametricEstimator( - "g_estimation_post_widder", {"x_data": x_data, "mu": mu, "sigma": sigma, "n": degree} + estimator = NMVSemiparEstim( + "density_estim_post_widder", {"x_data": x_data, "mu": mu, "sigma": sigma, "n": degree} ) est = estimator.estimate(sample) est_data = est.list_value diff --git a/tests/algorithms/nmv_algorithms/semiparametric_mu_estimation/__init__.py b/tests/procedures/nmv_procedures/semiparametric_mu_estimation/__init__.py similarity index 100% rename from tests/algorithms/nmv_algorithms/semiparametric_mu_estimation/__init__.py rename to tests/procedures/nmv_procedures/semiparametric_mu_estimation/__init__.py diff --git a/tests/algorithms/nmv_algorithms/semiparametric_mu_estimation/test_semiparametric_mu_estimation.py b/tests/procedures/nmv_procedures/semiparametric_mu_estimation/test_semiparametric_mu_estimation.py similarity index 88% rename from tests/algorithms/nmv_algorithms/semiparametric_mu_estimation/test_semiparametric_mu_estimation.py rename to tests/procedures/nmv_procedures/semiparametric_mu_estimation/test_semiparametric_mu_estimation.py index f1736f8..bcb9154 100644 --- a/tests/algorithms/nmv_algorithms/semiparametric_mu_estimation/test_semiparametric_mu_estimation.py +++ b/tests/procedures/nmv_procedures/semiparametric_mu_estimation/test_semiparametric_mu_estimation.py @@ -3,7 +3,7 @@ import pytest from scipy.stats import expon, gamma, halfnorm, pareto -from src.estimators.semiparametric.nmv_semiparametric_estimator import NMVSemiParametricEstimator +from src.estimators.semiparametric.nmv_semiparametric_estimator import NMVSemiparEstim from src.generators.nmv_generator import NMVGenerator from src.mixtures.nmv_mixture import NormalMeanVarianceMixtures @@ -15,7 +15,7 @@ class TestSemiParametricMuEstimation: def test_mu_estimation_expon_no_parameters(self, real_mu: float) -> None: mixture = NormalMeanVarianceMixtures("canonical", alpha=0, mu=real_mu, distribution=expon) sample = self.generator.canonical_generate(mixture, 10000) - estimator = NMVSemiParametricEstimator("mu_estimation") + estimator = NMVSemiparEstim("mu_estimation") est = estimator.estimate(sample) assert abs(real_mu - est.value) < 1 and est.success is True @@ -23,7 +23,7 @@ def test_mu_estimation_expon_no_parameters(self, real_mu: float) -> None: def test_mu_estimation_expon_no_parameters_small(self, real_mu: float) -> None: mixture = NormalMeanVarianceMixtures("canonical", alpha=0, mu=real_mu, distribution=expon) sample = self.generator.canonical_generate(mixture, 10000) - estimator = NMVSemiParametricEstimator("mu_estimation") + estimator = NMVSemiparEstim("mu_estimation") est = estimator.estimate(sample) assert abs(real_mu - est.value) < 1 and est.success is True @@ -31,7 +31,7 @@ def test_mu_estimation_expon_no_parameters_small(self, real_mu: float) -> None: def test_mu_estimation_pareto_no_parameters(self, real_mu: float) -> None: mixture = NormalMeanVarianceMixtures("canonical", alpha=0, mu=real_mu, distribution=pareto(2.62)) sample = self.generator.canonical_generate(mixture, 50000) - estimator = NMVSemiParametricEstimator("mu_estimation") + estimator = NMVSemiparEstim("mu_estimation") est = estimator.estimate(sample) assert abs(real_mu - est.value) < 1 and est.success is True @@ -39,7 +39,7 @@ def test_mu_estimation_pareto_no_parameters(self, real_mu: float) -> None: def test_mu_estimation_pareto_no_parameters_small(self, real_mu: float) -> None: mixture = NormalMeanVarianceMixtures("canonical", alpha=0, mu=real_mu, distribution=pareto(2.62)) sample = self.generator.canonical_generate(mixture, 50000) - estimator = NMVSemiParametricEstimator("mu_estimation") + estimator = NMVSemiparEstim("mu_estimation") est = estimator.estimate(sample) assert abs(real_mu - est.value) < 1 and est.success is True @@ -47,7 +47,7 @@ def test_mu_estimation_pareto_no_parameters_small(self, real_mu: float) -> None: def test_mu_estimation_halfnorm_no_parameters(self, real_mu: float) -> None: mixture = NormalMeanVarianceMixtures("canonical", alpha=0, mu=real_mu, distribution=halfnorm) sample = self.generator.canonical_generate(mixture, 10000) - estimator = NMVSemiParametricEstimator("mu_estimation") + estimator = NMVSemiparEstim("mu_estimation") est = estimator.estimate(sample) assert abs(real_mu - est.value) < 1 and est.success is True @@ -55,7 +55,7 @@ def test_mu_estimation_halfnorm_no_parameters(self, real_mu: float) -> None: def test_mu_estimation_gamma_no_parameters(self, real_mu: float) -> None: mixture = NormalMeanVarianceMixtures("canonical", alpha=0, mu=real_mu, distribution=gamma(2)) sample = self.generator.canonical_generate(mixture, 10000) - estimator = NMVSemiParametricEstimator("mu_estimation") + estimator = NMVSemiparEstim("mu_estimation") est = estimator.estimate(sample) assert abs(real_mu - est.value) < 1 and est.success is True @@ -64,7 +64,7 @@ def test_mu_estimation_expon_1_parameter_m_positive(self, params: dict) -> None: real_mu = 1 mixture = NormalMeanVarianceMixtures("canonical", alpha=0, mu=real_mu, distribution=expon) sample = self.generator.canonical_generate(mixture, 10000) - estimator = NMVSemiParametricEstimator("mu_estimation", params) + estimator = NMVSemiparEstim("mu_estimation", params) est = estimator.estimate(sample) assert abs(real_mu - est.value) < 1 and est.success is True @@ -73,7 +73,7 @@ def test_mu_estimation_expon_1_parameter_m_negative(self, params: dict) -> None: real_mu = -1 mixture = NormalMeanVarianceMixtures("canonical", alpha=0, mu=real_mu, distribution=expon) sample = self.generator.canonical_generate(mixture, 10000) - estimator = NMVSemiParametricEstimator("mu_estimation", params) + estimator = NMVSemiparEstim("mu_estimation", params) est = estimator.estimate(sample) assert abs(real_mu - est.value) < 1 and est.success is True @@ -84,7 +84,7 @@ def test_mu_estimation_expon_1_parameter_max_iterations(self, params: dict) -> N real_mu = 1 mixture = NormalMeanVarianceMixtures("canonical", alpha=0, mu=real_mu, distribution=expon) sample = self.generator.canonical_generate(mixture, 10000) - estimator = NMVSemiParametricEstimator("mu_estimation", params) + estimator = NMVSemiparEstim("mu_estimation", params) est = estimator.estimate(sample) assert abs(real_mu - est.value) < 1 and est.success is True @@ -93,7 +93,7 @@ def test_mu_estimation_expon_1_parameter_m_is_best_estimation(self, params: dict real_mu = 10 mixture = NormalMeanVarianceMixtures("canonical", alpha=0, mu=real_mu, distribution=expon) sample = self.generator.canonical_generate(mixture, 10000) - estimator = NMVSemiParametricEstimator("mu_estimation", params) + estimator = NMVSemiparEstim("mu_estimation", params) est = estimator.estimate(sample) assert abs(est.value == params["m"]) and est.success is False @@ -104,7 +104,7 @@ def test_mu_estimation_expon_2_parameters_tol_positive(self, params: dict) -> No real_mu = 1 mixture = NormalMeanVarianceMixtures("canonical", alpha=0, mu=real_mu, distribution=expon) sample = self.generator.canonical_generate(mixture, 10000) - estimator = NMVSemiParametricEstimator("mu_estimation", params) + estimator = NMVSemiparEstim("mu_estimation", params) est = estimator.estimate(sample) assert abs(real_mu - est.value) < 1 and est.success is True @@ -115,7 +115,7 @@ def test_mu_estimation_expon_2_parameters_tol_negative(self, params: dict) -> No real_mu = -1 mixture = NormalMeanVarianceMixtures("canonical", alpha=0, mu=real_mu, distribution=expon) sample = self.generator.canonical_generate(mixture, 10000) - estimator = NMVSemiParametricEstimator("mu_estimation", params) + estimator = NMVSemiparEstim("mu_estimation", params) est = estimator.estimate(sample) assert abs(real_mu - est.value) < 1 and est.success is True @@ -132,7 +132,7 @@ def test_mu_estimation_expon_3_parameters_omega_positive(self, params: dict) -> real_mu = 1 mixture = NormalMeanVarianceMixtures("canonical", alpha=0, mu=real_mu, distribution=expon) sample = self.generator.canonical_generate(mixture, 10000) - estimator = NMVSemiParametricEstimator("mu_estimation", params) + estimator = NMVSemiparEstim("mu_estimation", params) est = estimator.estimate(sample) assert abs(real_mu - est.value) < 1 and est.success is True @@ -149,7 +149,7 @@ def test_mu_estimation_expon_3_parameters_omega_negative(self, params: dict) -> real_mu = -1 mixture = NormalMeanVarianceMixtures("canonical", alpha=0, mu=real_mu, distribution=expon) sample = self.generator.canonical_generate(mixture, 10000) - estimator = NMVSemiParametricEstimator("mu_estimation", params) + estimator = NMVSemiparEstim("mu_estimation", params) est = estimator.estimate(sample) assert abs(real_mu - est.value) < 1 and est.success is True @@ -168,7 +168,7 @@ def test_mu_estimation_expon_3_parameters_all_positive(self, params: dict) -> No real_mu = 1 mixture = NormalMeanVarianceMixtures("canonical", alpha=0, mu=real_mu, distribution=expon) sample = self.generator.canonical_generate(mixture, 10000) - estimator = NMVSemiParametricEstimator("mu_estimation", params) + estimator = NMVSemiparEstim("mu_estimation", params) est = estimator.estimate(sample) assert abs(real_mu - est.value) < 1 and est.success is True @@ -187,6 +187,6 @@ def test_mu_estimation_expon_3_parameters_all_negative(self, params: dict) -> No real_mu = -1 mixture = NormalMeanVarianceMixtures("canonical", alpha=0, mu=real_mu, distribution=expon) sample = self.generator.canonical_generate(mixture, 10000) - estimator = NMVSemiParametricEstimator("mu_estimation", params) + estimator = NMVSemiparEstim("mu_estimation", params) est = estimator.estimate(sample) assert abs(real_mu - est.value) < 1 and est.success is True diff --git a/tests/algorithms/nmv_algorithms/semiparametric_mu_estimation/test_validate_kwargs.py b/tests/procedures/nmv_procedures/semiparametric_mu_estimation/test_validate_kwargs.py similarity index 82% rename from tests/algorithms/nmv_algorithms/semiparametric_mu_estimation/test_validate_kwargs.py rename to tests/procedures/nmv_procedures/semiparametric_mu_estimation/test_validate_kwargs.py index 781caf5..c0dc5de 100644 --- a/tests/algorithms/nmv_algorithms/semiparametric_mu_estimation/test_validate_kwargs.py +++ b/tests/procedures/nmv_procedures/semiparametric_mu_estimation/test_validate_kwargs.py @@ -3,7 +3,7 @@ import numpy as np import pytest -from src.algorithms.semiparam_algorithms.nvm_semi_param_algorithms.mu_estimation import SemiParametricMuEstimation +from src.procedures.semiparametric.nvm_semiparametric.mu_estimation import NMVEstimationMu def _test_omega(x: float) -> float: @@ -18,7 +18,7 @@ class TestValidateKwargs: ) def test_set_default_params_len_1_value_error_m(self, params: dict) -> None: with pytest.raises(ValueError, match="Expected positive integer as parameter m"): - SemiParametricMuEstimation()._validate_kwargs(**params) + NMVEstimationMu()._validate_kwargs(**params) @pytest.mark.parametrize( "params", @@ -26,12 +26,12 @@ def test_set_default_params_len_1_value_error_m(self, params: dict) -> None: ) def test_set_default_params_len_1_value_error_tolerance(self, params: dict) -> None: with pytest.raises(ValueError, match="Expected positive float as parameter tolerance"): - SemiParametricMuEstimation()._validate_kwargs(**params) + NMVEstimationMu()._validate_kwargs(**params) @pytest.mark.parametrize("params", [{"omega": 1}, {"omega": []}, {"omega": "str"}, {"omega": ()}]) def test_set_default_params_len_3_value_error_omega(self, params: dict) -> None: with pytest.raises(ValueError, match="Expected callable object as parameter omega"): - SemiParametricMuEstimation()._validate_kwargs(**params) + NMVEstimationMu()._validate_kwargs(**params) @pytest.mark.parametrize( "params", @@ -47,7 +47,7 @@ def test_set_default_params_len_3_value_error_omega(self, params: dict) -> None: ) def test_set_default_params_len_1_value_error_max_iterations(self, params: dict) -> None: with pytest.raises(ValueError, match="Expected positive integer as parameter max_iterations"): - SemiParametricMuEstimation()._validate_kwargs(**params) + NMVEstimationMu()._validate_kwargs(**params) @pytest.mark.parametrize( "params", @@ -62,7 +62,7 @@ def test_set_default_params_len_1_value_error_max_iterations(self, params: dict) ) def test_set_default_params_len_2_value_error_m(self, params: dict) -> None: with pytest.raises(ValueError, match="Expected positive integer as parameter m"): - SemiParametricMuEstimation()._validate_kwargs(**params) + NMVEstimationMu()._validate_kwargs(**params) @pytest.mark.parametrize( "params", @@ -73,7 +73,7 @@ def test_set_default_params_len_2_value_error_m(self, params: dict) -> None: ], ) def test_set_default_params_len_3_correct(self, params: dict) -> None: - SemiParametricMuEstimation()._validate_kwargs(**params) + NMVEstimationMu()._validate_kwargs(**params) @pytest.mark.parametrize( "params", @@ -84,17 +84,17 @@ def test_set_default_params_len_3_correct(self, params: dict) -> None: ], ) def test_set_default_params_len_4_correct(self, params: dict) -> None: - SemiParametricMuEstimation()._validate_kwargs(**params) + NMVEstimationMu()._validate_kwargs(**params) @pytest.mark.parametrize( "params", [{"m": 100, "tolerance": 1 / 10}, {"m": 1000, "tolerance": 10**-9}, {"m": 1, "tolerance": 1}] ) def test_set_default_params_len_2_correct(self, params: dict) -> None: - SemiParametricMuEstimation()._validate_kwargs(**params) + NMVEstimationMu()._validate_kwargs(**params) @pytest.mark.parametrize("params", [{"m": 100}, {"m": 1000}, {"m": 1}]) def test_set_default_params_len_1_correct(self, params: dict) -> None: - SemiParametricMuEstimation()._validate_kwargs(**params) + NMVEstimationMu()._validate_kwargs(**params) @pytest.mark.parametrize( "params", @@ -105,14 +105,14 @@ def test_set_default_params_len_1_correct(self, params: dict) -> None: ], ) def test_init_set_default_params_len_3_correct(self, params: dict) -> None: - SemiParametricMuEstimation(np.array([1]), **params) + NMVEstimationMu(np.array([1]), **params) @pytest.mark.parametrize( "params", [{"m": 100, "tolerance": 1 / 10}, {"m": 1000, "tolerance": 10**-9}, {"m": 1, "tolerance": 1}] ) def test_init_set_default_params_len_2_correct(self, params: dict) -> None: - SemiParametricMuEstimation(np.array([1]), **params) + NMVEstimationMu(np.array([1]), **params) @pytest.mark.parametrize("params", [{"m": 100}, {"m": 1000}, {"m": 1}]) def test_init_set_default_params_len_1_correct(self, params: dict) -> None: - SemiParametricMuEstimation(np.array([1]), **params) + NMVEstimationMu(np.array([1]), **params) \ No newline at end of file diff --git a/tests/algorithms/nv_algorithms/__init__.py b/tests/procedures/nv_procedures/__init__.py similarity index 100% rename from tests/algorithms/nv_algorithms/__init__.py rename to tests/procedures/nv_procedures/__init__.py diff --git a/tests/algorithms/nv_algorithms/test_g_estimation.py b/tests/procedures/nv_procedures/test_g_estimation.py similarity index 84% rename from tests/algorithms/nv_algorithms/test_g_estimation.py rename to tests/procedures/nv_procedures/test_g_estimation.py index f618be5..6371bb4 100644 --- a/tests/algorithms/nv_algorithms/test_g_estimation.py +++ b/tests/procedures/nv_procedures/test_g_estimation.py @@ -4,7 +4,7 @@ import pytest from scipy.stats import expon -from src.estimators.semiparametric.nv_semiparametric_estimator import NVSemiParametricEstimator +from src.estimators.semiparametric.nv_semiparametric_estimator import NVSemiparEstim from src.generators.nv_generator import NVGenerator from src.mixtures.nv_mixture import NormalVarianceMixtures @@ -17,8 +17,8 @@ def test_g_estimation_expon(self, x_data) -> None: n = 100 mixture = NormalVarianceMixtures("canonical", alpha=0, distribution=expon) sample = NVGenerator().canonical_generate(mixture, n) - estimator = NVSemiParametricEstimator( - "g_estimation_given_mu", {"x_data": x_data, "u_value": 7.6, "v_value": 0.9} + estimator = NVSemiparEstim( + "density_estim_inv_mellin_quad_int", {"x_data": x_data, "u_value": 7.6, "v_value": 0.9} ) est = estimator.estimate(sample) error = 0.0 diff --git a/tests/algorithms/nv_algorithms/test_post_widder_nv.py b/tests/procedures/nv_procedures/test_post_widder_nv.py similarity index 86% rename from tests/algorithms/nv_algorithms/test_post_widder_nv.py rename to tests/procedures/nv_procedures/test_post_widder_nv.py index e2fdee3..6dd917d 100644 --- a/tests/algorithms/nv_algorithms/test_post_widder_nv.py +++ b/tests/procedures/nv_procedures/test_post_widder_nv.py @@ -3,7 +3,7 @@ from mpmath import ln from scipy.stats import expon, gamma -from src.estimators.semiparametric.nv_semiparametric_estimator import NVSemiParametricEstimator +from src.estimators.semiparametric.nv_semiparametric_estimator import NVSemiparEstim from src.generators.nv_generator import NVGenerator from src.mixtures.nv_mixture import NormalVarianceMixtures @@ -20,8 +20,8 @@ def test_post_widder_expon(self, sigma, degree, sample_size) -> None: sample = NVGenerator().classical_generate(mixture, sample_size) x_data = np.linspace(0.5, 10.0, 30) - estimator = NVSemiParametricEstimator( - "g_estimation_post_widder", {"x_data": x_data, "sigma": sigma, "n": degree} + estimator = NVSemiparEstim( + "density_estim_post_widder", {"x_data": x_data, "sigma": sigma, "n": degree} ) est = estimator.estimate(sample) est_data = est.list_value @@ -43,8 +43,8 @@ def test_post_widder_gamma(self, sigma, degree, sample_size, a) -> None: sample = NVGenerator().classical_generate(mixture, sample_size) x_data = np.linspace(0.5, 10.0, 30) - estimator = NVSemiParametricEstimator( - "g_estimation_post_widder", {"x_data": x_data, "sigma": sigma, "n": degree} + estimator = NVSemiparEstim( + "density_estim_post_widder", {"x_data": x_data, "sigma": sigma, "n": degree} ) est = estimator.estimate(sample) est_data = est.list_value diff --git a/tests/algorithms/support_algorithms/test_rqmc.py b/tests/procedures/support_procedures/test_rqmc.py similarity index 97% rename from tests/algorithms/support_algorithms/test_rqmc.py rename to tests/procedures/support_procedures/test_rqmc.py index b5b6558..97171c1 100644 --- a/tests/algorithms/support_algorithms/test_rqmc.py +++ b/tests/procedures/support_procedures/test_rqmc.py @@ -3,7 +3,7 @@ import numpy as np import pytest -from src.algorithms.support_algorithms.rqmc import RQMC +from src.procedures.support.rqmc import RQMC def loss_func(true_func: Callable, rqms: Callable, count: int):